Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_ivp_2nd_rkn (d02la)

## Purpose

nag_ode_ivp_2nd_rkn (d02la) is a function for integrating a non-stiff system of second-order ordinary differential equations using Runge–Kutta–Nystrom techniques.

## Syntax

[t, y, yp, ydp, rwork, ifail] = d02la(fcn, t, tend, y, yp, ydp, rwork, 'neq', neq)
[t, y, yp, ydp, rwork, ifail] = nag_ode_ivp_2nd_rkn(fcn, t, tend, y, yp, ydp, rwork, 'neq', neq)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lrwork has been removed from the interface
.

## Description

Given the initial values x,y1,y2,,yneq,y1,y2,,yneq$x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}},{y}_{1}^{\prime },{y}_{2}^{\prime },\dots ,{y}_{{\mathbf{neq}}}^{\prime }$ nag_ode_ivp_2nd_rkn (d02la) integrates a non-stiff system of second-order differential equations of the type
 yi ′ ′ = fi(x,y1,y2, … ,yneq),  i = 1,2, … ,neq, $yi′′=fi(x,y1,y2,…,yneq), i=1,2,…,neq,$
from x = t$x={\mathbf{t}}$ to x = tend$x={\mathbf{tend}}$ using a Runge–Kutta–Nystrom formula pair. The system is defined by fcn, which evaluates fi${f}_{i}$ in terms of x$x$ and y1,y2,,yneq${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$, where y1,y2,,yneq${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$ are supplied at x$x$.
There are two Runge–Kutta–Nystrom formula pairs implemented in this function. The lower order method is intended if you have moderate accuracy requirements and may be used in conjunction with the interpolation function nag_ode_ivp_2nd_rkn_interp (d02lz) to produce solutions and derivatives at user-specified points. The higher order method is intended if you have high accuracy requirements.
In one-step mode the function returns approximations to the solution, derivative and fi${f}_{i}$ at each integration point. In interval mode these values are returned at the end of the integration range. You select the order of the method, the mode of operation, the error control and various optional inputs by a prior call to nag_ode_ivp_2nd_rkn_setup (d02lx).
For a description of the Runge–Kutta–Nystrom formula pairs see Dormand et al. (1986a) and Dormand et al. (1986b) and for a description of their practical implementation see Brankin et al. (1989).

## References

Brankin R W, Dormand J R, Gladwell I, Prince P J and Seward W L (1989) Algorithm 670: A Runge–Kutta–Nystrom Code ACM Trans. Math. Software 15 31–40
Dormand J R, El–Mikkawy M E A and Prince P J (1986a) Families of Runge–Kutta–Nystrom formulae Mathematical Report TPMR 86-1 Teesside Polytechnic
Dormand J R, El–Mikkawy M E A and Prince P J (1986b) High order embedded Runge–Kutta–Nystrom formulae Mathematical Report TPMR 86-2 Teesside Polytechnic

## Parameters

### Compulsory Input Parameters

1:     fcn – function handle or string containing name of m-file
fcn must evaluate the functions fi${f}_{i}$ (that is the second derivatives yi${y}_{i}^{\prime \prime }$) for given values of its arguments x$x$, y1,y2,,yneq${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$.
[f] = fcn(neq, t, y)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential equations.
2:     t – double scalar
x$x$, the value of the argument.
3:     y(neq) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the value of the argument.

Output Parameters

1:     f(neq) – double array
The value of fi${f}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
2:     t – double scalar
The initial value of the independent variable x$x$.
Constraint: ${\mathbf{t}}\ne {\mathbf{tend}}$.
3:     tend – double scalar
The end point of the range of integration. If ${\mathbf{tend}}<{\mathbf{t}}$ on initial entry, integration will proceed in the negative direction. tend may be reset, in the direction of integration, before any continuation call.
4:     y(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
The initial values of the solution y1,y2,,yneq${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$.
5:     yp(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
The initial values of the derivatives y1,y2,,yneq${y}_{1}^{\prime },{y}_{2}^{\prime },\dots ,{y}_{{\mathbf{neq}}}^{\prime }$.
6:     ydp(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
Must be unchanged from a previous call to nag_ode_ivp_2nd_rkn (d02la).
7:     rwork(lrwork) – double array
This must be the same parameter rwork as supplied to nag_ode_ivp_2nd_rkn_setup (d02lx). It is used to pass information from nag_ode_ivp_2nd_rkn_setup (d02lx) to nag_ode_ivp_2nd_rkn (d02la), and from nag_ode_ivp_2nd_rkn (d02la) to both nag_ode_ivp_2nd_rkn_diag (d02ly) and nag_ode_ivp_2nd_rkn_interp (d02lz). Therefore the contents of this array must not be changed before the call to nag_ode_ivp_2nd_rkn (d02la) or calling either of the functions nag_ode_ivp_2nd_rkn_diag (d02ly) and nag_ode_ivp_2nd_rkn_interp (d02lz).

### Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays y, yp, ydp. (An error is raised if these dimensions are not equal.)
The number of second-order ordinary differential equations to be solved by nag_ode_ivp_2nd_rkn (d02la). It must contain the same value as the parameter neq used in a prior call to nag_ode_ivp_2nd_rkn_setup (d02lx).
Constraint: neq1${\mathbf{neq}}\ge 1$.

lrwork

### Output Parameters

1:     t – double scalar
The value of the independent variable, which is usually tend, unless an error has occurred or the code is operating in one-step mode. If the integration is to be continued, possibly with a new value for tend, t must not be changed.
2:     y(neq) – double array
The computed values of the solution at the exit value of t. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.
3:     yp(neq) – double array
The computed values of the derivatives at the exit value of t. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.
4:     ydp(neq) – double array
The computed values of the second derivative at the exit value of t, unless illegal input is detected, in which case the elements of ydp may not have been initialized. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.
5:     rwork(lrwork) – double array
6:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
Illegal input detected, i.e., one of the following conditions:
• on any call, ${\mathbf{t}}={\mathbf{tend}}$, or the value of neq or lrwork has been altered;
• on a continuation call, the direction of integration has been changed;
• nag_ode_ivp_2nd_rkn_setup (d02lx) had not been called previously, or the previous call to nag_ode_ivp_2nd_rkn_setup (d02lx) resulted in an error exit.
This error exit can be caused if elements of rwork have been overwritten.
ifail = 2${\mathbf{ifail}}=2$
The maximum number of steps has been attempted. (See parameter maxstp in nag_ode_ivp_2nd_rkn_setup (d02lx).) If integration is to be continued then you need only reset ifail and call the function again and a further maxstp steps will be attempted.
ifail = 3${\mathbf{ifail}}=3$
In order to satisfy the error requirements, the step size needed is too small for the machine precision being used.
ifail = 4${\mathbf{ifail}}=4$
The code has detected two successive error exits at the current value of x$x$ and cannot proceed. Check all input variables.
ifail = 5${\mathbf{ifail}}=5$
The code has detected inefficient use of the integration method. The step size has been reduced by a significant amount too often in order to hit the output points specified by tend. (Of the last 100$100$ or more successful steps more than 10%$10%$ are steps with sizes that have had to be reduced by a factor of greater than a half.)

## Accuracy

The accuracy of integration is determined by the parameters tol, thres and thresp in a prior call to nag_ode_ivp_2nd_rkn_setup (d02lx). Note that only the local error at each step is controlled by these parameters. The error estimates obtained are not strict bounds but are usually reliable over one step. Over a number of steps the overall error may accumulate in various ways, depending on the system. The code is designed so that a reduction in tol should lead to an approximately proportional reduction in the error. You are strongly recommended to call nag_ode_ivp_2nd_rkn (d02la) with more than one value for tol and to compare the results obtained to estimate their accuracy.
The accuracy obtained depends on the type of error test used. If the solution oscillates around zero a relative error test should be avoided, whereas if the solution is exponentially increasing an absolute error test should not be used. For a description of the error test see the specifications of the parameters tol, thres and thresp in function document nag_ode_ivp_2nd_rkn_setup (d02lx).

If nag_ode_ivp_2nd_rkn (d02la) fails with ${\mathbf{ifail}}={\mathbf{3}}$ then the value of tol may be so small that a solution cannot be obtained, in which case the function should be called again with a larger value for tol. If the accuracy requested is really needed then you should consider whether there is a more fundamental difficulty. For example:
 (a) in the region of a singularity the solution components will usually be of a large magnitude. nag_ode_ivp_2nd_rkn (d02la) could be used in one-step mode to monitor the size of the solution with the aim of trapping the solution before the singularity. In any case numerical integration cannot be continued through a singularity, and analytical treatment may be necessary; (b) if the solution contains fast oscillatory components, the function will require a very small step size to preserve stability. This will usually be exhibited by excessive computing time and sometimes an error exit with ${\mathbf{ifail}}={\mathbf{3}}$. The Runge–Kutta–Nystrom methods are not efficient in such cases and you should consider reposing your problem as a system of first-order ordinary differential equations and then using a function from sub-chapter D02M–N with the Blend formulae (see nag_ode_ivp_stiff_dassl (d02mv)).
nag_ode_ivp_2nd_rkn (d02la) can be used for producing results at short intervals (for example, for tabulation), in two ways. By far the less efficient is to call nag_ode_ivp_2nd_rkn (d02la) successively over short intervals, t + (i1) × h$t+\left(i-1\right)×h$ to t + i × h$t+i×h$, although this is the only way if the higher order method has been selected and precisely not what it is intended for. A more efficient way, only for use when the lower order method has been selected, is to use nag_ode_ivp_2nd_rkn (d02la) in one-step mode. The output values of parameters y, yp, ydp, t and rwork are set correctly for a call to nag_ode_ivp_2nd_rkn_interp (d02lz) to compute the solution and derivative at the required points.

## Example

This example solves the following system (the two body problem)
 y1 ′ ′ = − y1 / (y12 + y22)3 / 2 y2 ′ ′ = − y2 / (y12 + y22)3 / 2
$y1′′ = -y1/(y12+y22)3/2 y2′′ = -y2/(y12+y22)3/2$
over the range [0,20]$\left[0,20\right]$ with initial conditions y1 = 1.0ε${y}_{1}=1.0-\epsilon$, y2 = 0.0${y}_{2}=0.0$, y1 = 0.0${y}_{1}^{\prime }=0.0$ and y2 = sqrt( ((1 + ε)/(1ε)) )${y}_{2}^{\prime }=\sqrt{\left(\frac{1+\epsilon }{1-\epsilon }\right)}$ where ε$\epsilon$, the eccentricity, is 0.5$0.5$. The system is solved using the lower order method with relative local error tolerances 1.0e−4$\text{1.0e−4}$ and 1.0e−5$\text{1.0e−5}$ and default threshold tolerances. nag_ode_ivp_2nd_rkn (d02la) is used in one-step mode (onestp = true${\mathbf{onestp}}=\mathbf{true}$) and nag_ode_ivp_2nd_rkn_interp (d02lz) provides solution values at intervals of 2.0$2.0$.
```function nag_ode_ivp_2nd_rkn_example
% Initialize variables and arrays.
ecc = 0.5;
tend = 20;
tstart = 0;
tinc = 2;

neq = 2;
lrwork = 16+20*neq;
nwant = neq;
thres = zeros(neq, 1);
thresp = zeros(neq, 1);

fprintf('nag_ode_ivp_2nd_rkn example program results \n\n');

for itol = 4:5
tol = 10^(-itol);
thres(1) = 0.0;
thresp(1) = 0.0;
h = 0.0;
maxstp = 0;
start = true;
high = false;
onestp = true;
rwork = zeros(56,1);
t = 0;
y = [1-ecc; 0];
yp = [0; sqrt((1+ecc)/(1-ecc))];
ydp = zeros(neq, 1);

% nag_ode_ivp_2nd_rkn_setup is a setup routine to be called prior to nag_ode_ivp_2nd_rkn.
[start, rwork, ifail] = ...
nag_ode_ivp_2nd_rkn_setup(h, tol, thres, thresp, int64(maxstp),...
start, onestp, high, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_setup returned with ifail = %1d ',ifail);
end

fprintf('Calculation with tol = %1.1e \n\n',tol);
fprintf('   T         Y(1)          Y(2) ');
t = tstart;
tnext = t + tinc;
fprintf('\n %4.1f ',t);
fprintf('   %10.5f',y);
fprintf('\n');

if itol == 4
ncall = 1;
tkeep = tnext;
ykeep = y';
end

while (t < tend)
[t, y, yp, ydp, rwork, ifail] = nag_ode_ivp_2nd_rkn(@fcn, t, tend, y,...
yp, ydp, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or integration problems.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn returned with ifail = %1d ',ifail);
end

while (tnext <= t)
% nag_ode_ivp_2nd_rkn_interp interpolates components of solution provided by nag_ode_ivp_2nd_rkn.
[ywant, ypwant, ifail] = ...
nag_ode_ivp_2nd_rkn_interp(t, y, yp, int64(neq), tnext,...
rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or interpolation problems.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_interp returned with ifail = %1d ',ifail);
end
ncall = ncall+1;

fprintf(' %4.1f ',tnext);
fprintf('   %10.5f',ywant);
fprintf('\n');

if itol == 4
tkeep(ncall) = tnext;
ykeep(ncall,:) = ywant';
end

tnext = tnext + tinc;
end
end

% nag_ode_ivp_2nd_rkn_diag is a diagnostic routine for use after calling nag_ode_ivp_2nd_rkn.
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] =...
nag_ode_ivp_2nd_rkn_diag(int64(neq), rwork);
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_diag returned with ifail = %1d ',ifail);
end

fprintf('\nNumber of successful steps = %1.0f\n',nsucc);
fprintf('Number of failed steps = %1.0f\n',nfail);
fprintf('\n\n');
end

% Now repeat the calculation with a smaller t increment, and store
% the trajectory in order to plot it.
tol = 10^(-4);
tinc = 0.1;
y = [1-ecc; 0];
yp = [0; sqrt((1+ecc)/(1-ecc))];

% nag_ode_ivp_2nd_rkn_setup is a setup routine to be called prior to nag_ode_ivp_2nd_rkn.
[start, rwork, ifail] = ...
nag_ode_ivp_2nd_rkn_setup(h, tol, thres, thresp, int64(maxstp),...
start, onestp, high, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_setup returned with ifail = %1d ',ifail);
end

t = tstart;
tnext = t + tinc;

nkeep = 1;
ykeep = y';

while (t < tend)
[t, y, yp, ydp, rwork, ifail] = nag_ode_ivp_2nd_rkn(@fcn, t, tend, y,...
yp, ydp, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or integration problems.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn returned with ifail = %1d ',ifail);
end

while (tnext <= t)
% nag_ode_ivp_2nd_rkn_interp interpolates components of solution provided by nag_ode_ivp_2nd_rkn.
[ywant, ypwant, ifail] = nag_ode_ivp_2nd_rkn_interp(t, y, yp, int64(neq), tnext,...
rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or interpolation problems.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_interp returned with ifail = %1d ',ifail);
end

nkeep = nkeep+1;
ykeep(nkeep,:) = ywant';

tnext = tnext + tinc;
end
end

% nag_ode_ivp_2nd_rkn_diag is a diagnostic routine for use after calling nag_ode_ivp_2nd_rkn.
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] =...
nag_ode_ivp_2nd_rkn_diag(int64(neq), rwork);
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: nag_ode_ivp_2nd_rkn_diag returned with ifail = %1d ',ifail);
end
% Find the index of the point halfway through the first orbit (i.e. the
% first point at which the y coord is negative).  This is used to calculate
% the end points of the orbits in display_plot.
nhalf = 1;
while ykeep(nhalf,2) >= 0
nhalf = nhalf+1;
end
% Plot results.
fig = figure('Number', 'off');
display_plot(ykeep, nhalf);

function ydp = fcn(neq, t, y)
% Evaluate second derivatives.
r = sqrt(y(1)^2+y(2)^2)^3;
f = zeros(2,1);
ydp(1) = -y(1)/r;
ydp(2) = -y(2)/r;
function display_plot(ykeep, nhalf)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.

% Find the end points of the three orbits.
norb1 = 2*(nhalf-1);
norb2 = 4*(nhalf-1);
norb3 = 6*(nhalf-1);

% Plot the orbits, including displacements.
plot(ykeep(1:norb1,1), ykeep(1:norb1,2), '-+', ...
ykeep(norb1+1:norb2,1), ykeep(norb1+1:norb2,2)+0.1, '--x', ...
ykeep(norb2+1:norb3,1), ykeep(norb2+1:norb3,2)+0.2, ':*');
title({'Second-order ODE Solution using Runge-Kutta-Nystrom.', ...
'The Two-body Problem (using shifts to distinguish orbits)'}, ...
titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('y', labFmt{:});
legend('1st orbit', '2nd orbit+(0,0.1)', '3rd orbit+(0,0.2)', ...
'Location', 'Best');
```
```
nag_ode_ivp_2nd_rkn example program results

Calculation with tol = 1.0e-04

T         Y(1)          Y(2)
0.0       0.50000      0.00000
2.0      -1.20573      0.61357
4.0      -1.33476     -0.47685
6.0       0.35748     -0.44558
8.0      -1.03762      0.73022
10.0      -1.42617     -0.32658
12.0       0.05515     -0.72032
14.0      -0.82880      0.81788
16.0      -1.48103     -0.16788
18.0      -0.26719     -0.84223
20.0      -0.57803      0.86339

Number of successful steps = 108
Number of failed steps = 16

Calculation with tol = 1.0e-05

T         Y(1)          Y(2)
0.0       0.50000      0.00000
2.0      -1.20573      0.61357
4.0      -1.33476     -0.47685
6.0       0.35748     -0.44558
8.0      -1.03762      0.73022
10.0      -1.42617     -0.32658
12.0       0.05516     -0.72031
14.0      -0.82880      0.81787
16.0      -1.48103     -0.16789
18.0      -0.26718     -0.84223
20.0      -0.57804      0.86338

Number of successful steps = 169
Number of failed steps = 15

```
```function d02la_example
% Initialize variables and arrays.
ecc = 0.5;
tend = 20;
tstart = 0;
tinc = 2;

neq = 2;
lrwork = 16+20*neq;
nwant = neq;
thres = zeros(neq, 1);
thresp = zeros(neq, 1);

fprintf('d02la example program results \n\n');

for itol = 4:5
tol = 10^(-itol);
thres(1) = 0.0;
thresp(1) = 0.0;
h = 0.0;
maxstp = 0;
start = true;
high = false;
onestp = true;
rwork = zeros(56,1);
t = 0;
y = [1-ecc; 0];
yp = [0; sqrt((1+ecc)/(1-ecc))];
ydp = zeros(neq, 1);

% d02lx is a setup routine to be called prior to d02la.
[start, rwork, ifail] = d02lx(h, tol, thres, thresp, int64(maxstp),...
start, onestp, high, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: d02lx returned with ifail = %1d ',ifail);
end

fprintf('Calculation with tol = %1.1e \n\n',tol);
fprintf('   T         Y(1)          Y(2) ');
t = tstart;
tnext = t + tinc;
fprintf('\n %4.1f ',t);
fprintf('   %10.5f',y);
fprintf('\n');

if itol == 4
ncall = 1;
tkeep = tnext;
ykeep = y';
end

while (t < tend)
[t, y, yp, ydp, rwork, ifail] = d02la(@fcn, t, tend, y,...
yp, ydp, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or integration problems.  Print message and exit.
error('Warning: d02la returned with ifail = %1d ',ifail);
end

while (tnext <= t)
% d02lz interpolates components of solution provided by d02la.
[ywant, ypwant, ifail] = d02lz(t, y, yp, int64(neq), tnext,...
rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or interpolation problems.  Print message and exit.
error('Warning: d02lz returned with ifail = %1d ',ifail);
end
ncall = ncall+1;

fprintf(' %4.1f ',tnext);
fprintf('   %10.5f',ywant);
fprintf('\n');

if itol == 4
tkeep(ncall) = tnext;
ykeep(ncall,:) = ywant';
end

tnext = tnext + tinc;
end
end

% d02ly is a diagnostic routine for use after calling d02la.
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] =...
d02ly(int64(neq), rwork);
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: d02ly returned with ifail = %1d ',ifail);
end

fprintf('\nNumber of successful steps = %1.0f\n',nsucc);
fprintf('Number of failed steps = %1.0f\n',nfail);
fprintf('\n\n');
end

% Now repeat the calculation with a smaller t increment, and store
% the trajectory in order to plot it.
tol = 10^(-4);
tinc = 0.1;
y = [1-ecc; 0];
yp = [0; sqrt((1+ecc)/(1-ecc))];

% d02lx is a setup routine to be called prior to d02la.
[start, rwork, ifail] = d02lx(h, tol, thres, thresp, int64(maxstp),...
start, onestp, high, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: d02lx returned with ifail = %1d ',ifail);
end

t = tstart;
tnext = t + tinc;

nkeep = 1;
ykeep = y';

while (t < tend)
[t, y, yp, ydp, rwork, ifail] = d02la(@fcn, t, tend, y,...
yp, ydp, rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or integration problems.  Print message and exit.
error('Warning: d02la returned with ifail = %1d ',ifail);
end

while (tnext <= t)
% d02lz interpolates components of solution provided by d02la.
[ywant, ypwant, ifail] = d02lz(t, y, yp, int64(neq), tnext,...
rwork, 'neq', int64(neq));
if ifail ~= 0
% Illegal input or interpolation problems.  Print message and exit.
error('Warning: d02lz returned with ifail = %1d ',ifail);
end

nkeep = nkeep+1;
ykeep(nkeep,:) = ywant';

tnext = tnext + tinc;
end
end

% d02ly is a diagnostic routine for use after calling d02la.
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] =...
d02ly(int64(neq), rwork);
if ifail ~= 0
% Illegal input.  Print message and exit.
error('Warning: d02ly returned with ifail = %1d ',ifail);
end
% Find the index of the point halfway through the first orbit (i.e. the
% first point at which the y coord is negative).  This is used to calculate
% the end points of the orbits in display_plot.
nhalf = 1;
while ykeep(nhalf,2) >= 0
nhalf = nhalf+1;
end
% Plot results.
fig = figure('Number', 'off');
display_plot(ykeep, nhalf);

function ydp = fcn(neq, t, y)
% Evaluate second derivatives.
r = sqrt(y(1)^2+y(2)^2)^3;
f = zeros(2,1);
ydp(1) = -y(1)/r;
ydp(2) = -y(2)/r;
function display_plot(ykeep, nhalf)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.

% Find the end points of the three orbits.
norb1 = 2*(nhalf-1);
norb2 = 4*(nhalf-1);
norb3 = 6*(nhalf-1);

% Plot the orbits, including displacements.
plot(ykeep(1:norb1,1), ykeep(1:norb1,2), '-+', ...
ykeep(norb1+1:norb2,1), ykeep(norb1+1:norb2,2)+0.1, '--x', ...
ykeep(norb2+1:norb3,1), ykeep(norb2+1:norb3,2)+0.2, ':*');
title({'Second-order ODE Solution using Runge-Kutta-Nystrom.', ...
'The Two-body Problem (using shifts to distinguish orbits)'}, ...
titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('y', labFmt{:});
legend('1st orbit', '2nd orbit+(0,0.1)', '3rd orbit+(0,0.2)', ...
'Location', 'Best');
```
```
d02la example program results

Calculation with tol = 1.0e-04

T         Y(1)          Y(2)
0.0       0.50000      0.00000
2.0      -1.20573      0.61357
4.0      -1.33476     -0.47685
6.0       0.35748     -0.44558
8.0      -1.03762      0.73022
10.0      -1.42617     -0.32658
12.0       0.05515     -0.72032
14.0      -0.82880      0.81788
16.0      -1.48103     -0.16788
18.0      -0.26719     -0.84223
20.0      -0.57803      0.86339

Number of successful steps = 108
Number of failed steps = 16

Calculation with tol = 1.0e-05

T         Y(1)          Y(2)
0.0       0.50000      0.00000
2.0      -1.20573      0.61357
4.0      -1.33476     -0.47685
6.0       0.35748     -0.44558
8.0      -1.03762      0.73022
10.0      -1.42617     -0.32658
12.0       0.05516     -0.72031
14.0      -0.82880      0.81787
16.0      -1.48103     -0.16789
18.0      -0.26718     -0.84223
20.0      -0.57804      0.86338

Number of successful steps = 169
Number of failed steps = 15

```