NAG Toolbox: nag_ode_ivp_stiff_exp_revcom (d02nm)

Purpose

Syntax

Description

% 
% data initialization
% 
call linear algebra setup routine 
call integrator setup routine 
irevcm = int32(0);
[..., irevcm, ...] = d02nm();
while (irevcm > 0)
  if (irevcm == 8)
    supply the jacobian matrix                                  (i)
  elseif (irevcm == 9)
    perform monitoring tasks requested by the user             (ii)
  elseif (irecvm == 1 or irevcm >= 3 and irevcm <= 5)
    evaluate the derivative                                   (iii)
  elseif (irevcm == 10)
    indicates an unsuccessful step 
  end 
  [..., irevcm, ...] = d02nm();
end
% 
% post processing (optional linear algebra diagnostic call 
% (sparse case only), optional integrator diagnostic call) 
%

[errloc, ifail] = d02za(rwork(51+neq:51+neq+neq-1), rwork(51:51+neq-1));
% Check ifail before proceeding

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{t}$ – double scalar

On initial entry: $t$, the value of the independent variable. The input value of t is used only on the first call as the initial point of the integration.

2:     $\mathrm{tout}$ – double scalar

On initial entry: the next value of $t$ at which a computed solution is desired. For the initial $t$, the input value of tout is used to determine the direction of integration. Integration is permitted in either direction (see also itask).

Constraint: $\mathbf{tout}\ne \mathbf{t}$.

3:     $\mathrm{y}\left(\mathbf{neq}\right)$ – double array

On initial entry: the values of the dependent variables (solution). On the first call the first neq elements of $y$ must contain the vector of initial values.

4:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

On intermediate re-entry: must be set to the derivatives as defined under the description of irevcm.

5:     $\mathrm{rwork}\left(50+4\times \mathbf{neq}\right)$ – double array

On initial entry: must be the same array as used by one of the method setup functions nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw), and by one of the storage setup functions nag_ode_ivp_stiff_bandjac_setup (d02nt), nag_ode_ivp_stiff_sparjac_setup (d02nu) or nag_ode_ivp_stiff_bdf (d02nv). The contents of rwork must not be changed between any call to a setup function and the first call to nag_ode_ivp_stiff_exp_revcom (d02nm).

On intermediate re-entry: elements of rwork must be set to quantities as defined under the description of irevcm.

6:     $\mathrm{rtol}\left(:\right)$ – double array

The dimension of the array rtol must be at least $1$ if $\mathbf{itol}=1$ or $2$, and at least $\mathbf{neq}$ otherwise

On initial entry: the relative local error tolerance.

Constraint: $\mathbf{rtol}\left(i\right)\ge 0.0$ for all relevant $i$ (see itol).

7:     $\mathrm{atol}\left(:\right)$ – double array

The dimension of the array atol must be at least $1$ if $\mathbf{itol}=1$ or $3$, and at least $\mathbf{neq}$ otherwise

On initial entry: the absolute local error tolerance.

Constraint: $\mathbf{atol}\left(i\right)\ge 0.0$ for all relevant $i$ (see itol).

8:     $\mathrm{itol}$ – int64int32nag_int scalar

On initial entry: a value to indicate the form of the local error test. itol indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) whether to interpret either or both of rtol or atol as a vector or a scalar. The error test to be satisfied is $‖e_i/w_i‖<1.0$, where $w_i$ is defined as follows:

itol	rtol	atol	$w_i$
1	scalar	scalar	$\mathbf{rtol}\left(1\right)\times \left|y_i\right|+\mathbf{atol}\left(1\right)$
2	scalar	vector	$\mathbf{rtol}\left(1\right)\times \left|y_i\right|+\mathbf{atol}\left(i\right)$
3	vector	scalar	$\mathbf{rtol}\left(i\right)\times \left|y_i\right|+\mathbf{atol}\left(1\right)$
4	vector	vector	$\mathbf{rtol}\left(i\right)\times \left|y_i\right|+\mathbf{atol}\left(i\right)$

$e_i$ is an estimate of the local error in $y_i$, computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup function.

Constraint: $\mathbf{itol}=1$, $2$, $3$ or $4$.

9:     $\mathrm{inform}\left(23\right)$ – int64int32nag_int array

10:   $\mathrm{ysav}\left(\mathit{ldysav},\mathbf{sdysav}\right)$ – double array

ldysav, the first dimension of the array, must satisfy the constraint $\mathit{ldysav}\ge \mathbf{neq}$.

On initial entry: the second dimension of the array ysav. an appropriate value for sdysav is described in the specifications of the integrator setup functions nag_ode_ivp_stiff_bdf (d02nv) and nag_ode_ivp_stiff_blend (d02nw). This value must be the same as that supplied to the integrator setup function.

11:   $\mathrm{wkjac}\left(\mathbf{nwkjac}\right)$ – double array

On intermediate re-entry: elements of the Jacobian as defined under the description of irevcm. If a numerical Jacobian was requested then wkjac is used for workspace.

12:   $\mathrm{jacpvt}\left(\mathit{njcpvt}\right)$ – int64int32nag_int array

On initial entry: the dimension of the array jacpvt. the actual size depends on the linear algebra method used. An appropriate value for njcpvt is described in the specifications of the linear algebra setup functions nag_ode_ivp_stiff_bandjac_setup (d02nt) and nag_ode_ivp_stiff_sparjac_setup (d02nu) for banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup function. When full matrix linear algebra is chosen, the array jacpvt is not used and hence njcpvt should be set to $1$.

13:   $\mathrm{imon}$ – int64int32nag_int scalar

On intermediate re-entry: may be reset to determine subsequent action in nag_ode_ivp_stiff_exp_revcom (d02nm).

$\mathbf{imon}=-2$

Integration is to be halted. A return will be made from nag_ode_ivp_stiff_exp_revcom (d02nm) to the calling (sub)program with $\mathbf{ifail}=\mathbf{12}$.

$\mathbf{imon}=-1$

Allow nag_ode_ivp_stiff_exp_revcom (d02nm) to continue with its own internal strategy. The integrator will try up to three restarts unless $\mathbf{imon}\ne -1$.

$\mathbf{imon}=0$

Return to the internal nonlinear equation solver, where the action taken is determined by the value of inln.

$\mathbf{imon}=1$

Normal exit to nag_ode_ivp_stiff_exp_revcom (d02nm) to continue integration.

$\mathbf{imon}=2$

Restart the integration at the current time point. The integrator will restart from order $1$ when this option is used. The solution y, provided by the monitr operation (see Description), will be used for the initial conditions.

$\mathbf{imon}=3$

Try to continue with the same step size and order as was to be used before entering the monitr operation (see Description). hmin and hmax may be altered if desired.

$\mathbf{imon}=4$

Continue the integration but using a new value of hnext and possibly new values of hmin and hmax.

14:   $\mathrm{inln}$ – int64int32nag_int scalar

On intermediate re-entry: with $\mathbf{imon}=0$ and $\mathbf{irevcm}=9$, inln specifies the action to be taken by the internal nonlinear equation solver. By setting $\mathbf{inln}=3$ and returning to nag_ode_ivp_stiff_exp_revcom (d02nm), the residual vector is evaluated and placed in $\mathbf{rwork}\left(50+2\times \mathbf{neq}+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$ and then the monitr operation (see Description) is invoked again. At present this is the only option available: inln must not be set to any other value.

15:   $\mathrm{ires}$ – int64int32nag_int scalar

On intermediate re-entry: should be unchanged unless one of the following actions is required of nag_ode_ivp_stiff_exp_revcom (d02nm) in which case ires should be set accordingly.

$\mathbf{ires}=2$

Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) that control should be passed back immediately to the calling (sub)program with the error indicator set to $\mathbf{ifail}=\mathbf{11}$.

$\mathbf{ires}=3$

Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) that an error condition has occurred in the solution vector, its time derivative or in the value of $t$. The integrator will use a smaller time step to try to avoid this condition. If this is not possible nag_ode_ivp_stiff_exp_revcom (d02nm) returns to the calling (sub)program with the error indicator set to $\mathbf{ifail}=\mathbf{7}$.

$\mathbf{ires}=4$

Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) to stop its current operation and to enter the monitr operation (see Description) immediately.

16:   $\mathrm{irevcm}$ – int64int32nag_int scalar

On initial entry: must contain $0$.

On intermediate re-entry: should remain unchanged.

Constraint: $\mathbf{irevcm}=0$, $1$, $3$, $4$, $5$, $8$, $9$ or $10$.

17:   $\mathrm{itask}$ – int64int32nag_int scalar

On initial entry: the task to be performed by the integrator.

$\mathbf{itask}=1$

Normal computation of output values of $y\left(t\right)$ at $t=\mathbf{tout}$ (by overshooting and interpolating).

$\mathbf{itask}=2$

Take one step only and return.

$\mathbf{itask}=3$

Stop at the first internal integration point at or beyond $t=\mathbf{tout}$ and return.

$\mathbf{itask}=4$

Normal computation of output values of $y\left(t\right)$ at $t=\mathbf{tout}$ but without overshooting $t=\mathbf{tcrit}$. tcrit must be specified as an option in one of the integrator setup functions before the first call to the integrator, or specified in the optional input function before a continuation call. tcrit (e.g., see nag_ode_ivp_stiff_bdf (d02nv)) may be equal to or beyond tout, but not before it in the direction of integration.

$\mathbf{itask}=5$

Take one step only and return, without passing tcrit (e.g., see nag_ode_ivp_stiff_bdf (d02nv)). tcrit must be specified under $\mathbf{itask}=4$.

Constraint: $1\le \mathbf{itask}\le 5$.

18:   $\mathrm{itrace}$ – int64int32nag_int scalar

On initial entry: the level of output that is printed by the integrator. itrace may take the value $-1$, $0$, $1$, $2$ or $3$.

$\mathbf{itrace}<-1$

$-1$ is assumed and similarly if $\mathbf{itrace}>3$, then $3$ is assumed.

$\mathbf{itrace}=-1$

No output is generated.

$\mathbf{itrace}=0$

Only warning messages are printed on the current error message unit (see nag_file_set_unit_error (x04aa)).

$\mathbf{itrace}>0$

Warning messages are printed as above, and on the current advisory message unit (see nag_file_set_unit_advisory (x04ab)) output is generated which details Jacobian entries, the nonlinear iteration and the time integration. The advisory messages are given in greater detail the larger the value of itrace.

Optional Input Parameters

1:     $\mathrm{neq}$ – int64int32nag_int scalar

Default: the dimension of the arrays y, ydot and the first dimension of the array ysav. (An error is raised if these dimensions are not equal.)

On initial entry: the number of differential equations to be solved.

Constraint: $\mathbf{neq}\ge 1$.

2:     $\mathrm{sdysav}$ – int64int32nag_int scalar

Default: the second dimension of the array ysav.

On initial entry: the second dimension of the array ysav. an appropriate value for sdysav is described in the specifications of the integrator setup functions nag_ode_ivp_stiff_bdf (d02nv) and nag_ode_ivp_stiff_blend (d02nw). This value must be the same as that supplied to the integrator setup function.

3:     $\mathrm{nwkjac}$ – int64int32nag_int scalar

Default: the dimension of the array wkjac.

On initial entry: the dimension of the array wkjac. the actual size depends on the linear algebra method used. An appropriate value for nwkjac is described in the specifications of the linear algebra setup functions nag_ode_ivp_stiff_fulljac_setup (d02ns), nag_ode_ivp_stiff_bandjac_setup (d02nt) and nag_ode_ivp_stiff_sparjac_setup (d02nu) for full, banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup function.

Output Parameters

1:     $\mathrm{t}$ – double scalar

On final exit: the value at which the computed solution $y$ is returned (usually at tout).

2:     $\mathrm{y}\left(\mathbf{neq}\right)$ – double array

On final exit: the computed solution vector evaluated at t (usually $t=\mathbf{tout}$).

3:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

On final exit: the time derivatives $y^{\prime }$ of the vector $y$ at the last integration point.

4:     $\mathrm{rwork}\left(50+4\times \mathbf{neq}\right)$ – double array

On intermediate exit: contains information for JAC, FCN and MONITR operations as described in Description and the argument irevcm.

5:     $\mathrm{inform}\left(23\right)$ – int64int32nag_int array

6:     $\mathrm{ysav}\left(\mathit{ldysav},\mathbf{sdysav}\right)$ – double array

7:     $\mathrm{wkjac}\left(\mathbf{nwkjac}\right)$ – double array

On intermediate exit: the Jacobian is overwritten.

8:     $\mathrm{jacpvt}\left(\mathit{njcpvt}\right)$ – int64int32nag_int array

9:     $\mathrm{imon}$ – int64int32nag_int scalar

On intermediate exit: used to pass information between nag_ode_ivp_stiff_exp_revcom (d02nm) and the MONITR operation (see Description). With $\mathbf{irevcm}=9$, imon contains a flag indicating under what circumstances the return from nag_ode_ivp_stiff_exp_revcom (d02nm) occurred:

$\mathbf{imon}=-2$

Exit from nag_ode_ivp_stiff_exp_revcom (d02nm) after $\mathbf{ires}=4$ caused an early termination (this facility could be used to locate discontinuities).

$\mathbf{imon}=-1$

The current step failed repeatedly.

$\mathbf{imon}=0$

Exit from nag_ode_ivp_stiff_exp_revcom (d02nm) after a call to the internal nonlinear equation solver.

$\mathbf{imon}=1$

The current step was successful.

10:   $\mathrm{inln}$ – int64int32nag_int scalar

On intermediate exit: contains a flag indicating the action to be taken, if any, by the internal nonlinear equation solver.

11:   $\mathrm{ires}$ – int64int32nag_int scalar

On intermediate exit: with $\mathbf{irevcm}=1$, $2$, $3$, $4$ or $5$, ires contains the value $1$.

12:   $\mathrm{irevcm}$ – int64int32nag_int scalar

On intermediate exit: indicates what action you must take before re-entering. The possible exit values of irevcm are $1$, $3$, $4$, $5$, $8$, $9$ or $10$, which should be interpreted as follows:

$\mathbf{irevcm}=1$, $3$, $4$ and $5$

Indicates that an FCN operation (see Description) is required: $y^{\prime }=g\left(t,y\right)$ must be supplied, where $\mathbf{y}\left(\mathit{i}\right)$ is located in $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$.

For $\mathbf{irevcm}=1$ or $3$, $y_{\mathit{i}}^{\prime }$ should be placed in location $\mathbf{rwork}\left(50+2\times \mathbf{neq}+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$.

For $\mathbf{irevcm}=4$, $y_{\mathit{i}}^{\prime }$ should be placed in location $\mathbf{rwork}\left(50+\mathbf{neq}+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$.

For $\mathbf{irevcm}=5$, $y_{\mathit{i}}^{\prime }$ should be placed in location $\mathbf{ydot}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$.

$\mathbf{irevcm}=8$

Indicates that a JAC operation (see Description) is required: the Jacobian matrix must be supplied.

If full matrix linear algebra is being used, then the $\left(i,j\right)$th element of the Jacobian must be stored in $\mathbf{wkjac}\left(\left(j-1\right)\times \mathbf{neq}+i\right)$.

If banded matrix linear algebra is being used then the $\left(i,j\right)$th element of the Jacobian
must be stored in $\mathbf{wkjac}\left(\left(i-1\right)\times m_B+k\right)$, where $m_B=m_L+m_U+1$ and
$k=\min \left(m_L-i+1,0\right)+j$; here $m_L$ and $m_U$ are the number of subdiagonals and superdiagonals, respectively, in the band.

If sparse matrix linear algebra is being used then nag_ode_ivp_stiff_sparjac_enq (d02nr) must be called to determine which column of the Jacobian is required and where it should be stored.

 [j, iplace] = d02nr(inform); 

will return in j the number of the column of the Jacobian that is required and will set $\mathbf{iplace}=1$ or $2$. If $\mathbf{iplace}=1$, then the $\left(i,j\right)$th element of the Jacobian must be stored in $\mathbf{rwork}\left(50+2\times \mathbf{neq}+i\right)$; otherwise it must be stored in $\mathbf{rwork}\left(50+\mathbf{neq}+i\right)$.

$\mathbf{irevcm}=9$

Indicates that a MONITR operation (see Description) can be performed.

$\mathbf{irevcm}=10$

Indicates that the current step was not successful, due to error test failure or convergence test failure. The only information supplied to you on this return is the current value of the independent variable $t$, located in $\mathbf{rwork}\left(19\right)$. No values must be changed before re-entering nag_ode_ivp_stiff_exp_revcom (d02nm); this facility enables you to determine the number of unsuccessful steps.

On final exit: $\mathbf{irevcm}=0$ indicated the user-specified task has been completed or an error has been encountered (see the descriptions for itask and ifail).

13:   $\mathrm{ifail}$ – int64int32nag_int scalar

On final exit: $\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

   $\mathbf{ifail}=1$

On entry, the integrator detected an illegal input, or that a linear algebra and/or integrator setup function has not been called prior to the call to the integrator. If $\mathbf{itrace}\ge 0$, the form of the error will be detailed on the current error message unit (see nag_file_set_unit_error (x04aa)).

   $\mathbf{ifail}=2$

The maximum number of steps specified has been taken (see the description of optional inputs in the integrator setup functions and the optional input continuation function, nag_ode_ivp_stiff_contin (d02nz)).

W  $\mathbf{ifail}=3$

With the given values of rtol and atol no further progress can be made across the integration range from the current point t. The components $\mathbf{y}\left(1\right),\mathbf{y}\left(2\right),\dots ,\mathbf{y}\left(\mathbf{neq}\right)$ contain the computed values of the solution at the current point t.

W  $\mathbf{ifail}=4$

There were repeated error test failures on an attempted step, before completing the requested task, but the integration was successful as far as t. The problem may have a singularity, or the local error requirements may be inappropriate.

W  $\mathbf{ifail}=5$

There were repeated convergence test failures on an attempted step, before completing the requested task, but the integration was successful as far as t. This may be caused by an inaccurate Jacobian matrix or one which is incorrectly computed.

W  $\mathbf{ifail}=6$

Some error weight $w_i$ became zero during the integration (see the description of itol). Pure relative error control $\left(\mathbf{atol}\left(i\right)=0.0\right)$ was requested on a variable (the $i$th) which has now vanished. The integration was successful as far as t.

   $\mathbf{ifail}=7$

The FCN operation (see Description) set the error flag $\mathbf{ires}=3$ continually despite repeated attempts by the integrator to avoid this.

   $\mathbf{ifail}=8$

Not used for this integrator.

   $\mathbf{ifail}=9$

A singular Jacobian $\frac{\partial r}{\partial y}$ has been encountered. This error exit is unlikely to be taken when solving explicit ordinary differential equations. You should check the problem formulation and Jacobian calculation.

   $\mathbf{ifail}=10$

An error occurred during Jacobian formulation or back-substitution (a more detailed error description may be directed to the current error message unit, see nag_file_set_unit_error (x04aa)).

W  $\mathbf{ifail}=11$

The FCN operation (see Description) signalled the integrator to halt the integration and return by setting $\mathbf{ires}=2$. Integration was successful as far as t.

W  $\mathbf{ifail}=12$

The MONITR operation (see Description) set $\mathbf{imon}=-2$ and so forced a return but the integration was successful as far as t.

W  $\mathbf{ifail}=13$

The requested task has been completed, but it is estimated that a small change in rtol and atol is unlikely to produce any change in the computed solution. (Only applies when you are not operating in one step mode, that is when $\mathbf{itask}\ne 2$ or $5$.)

   $\mathbf{ifail}=14$

The values of rtol and atol are so small that nag_ode_ivp_stiff_exp_revcom (d02nm) is unable to start the integration.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02nm_example

function d02nm_example


fprintf('d02nm example results\n\n');

% Initialize setup setup variables and arrays.
neq    = int64(3);
neqmax = neq;
maxord = int64(5);
sdysav = maxord+1;
petzld = false;
tcrit  = 0;
hmin   = 1.0e-10;
hmax   = 10;
h0     = 0;
maxstp = int64(200);
mxhnil = int64(5);
const  = zeros(6, 1);
rwork  = zeros(50+4*neq, 1);

% BDF with Newton iterations.
[const, rwork, ifail] = d02nv( ...
                               neqmax, sdysav, maxord, 'Newton', petzld, ...
                               const, tcrit, hmin, hmax, h0, maxstp, ...
                               mxhnil, 'Average-L2', rwork);
% Numerical Jacobian
nwkjac = neqmax*(neqmax + 1);
[rwork, ifail] = d02ns( ...
                        neq, neqmax, 'Numerical', nwkjac, rwork);

% Initialize variables and arrays.
njcpvt = int64(1);
wkjac  = zeros(nwkjac, 1);
ydot   = zeros(neq, 1);
ysave  = zeros(neq, sdysav);
inform(1:23) = int64(0);
jacpvt(1) = int64(0);
algequ = zeros(neq, 1);

% Integrate to tout by overshooting (itask = 1).
% At monitoring stages output intermediate solutions after interpolating.
t     = 0;
tout  = 10.0;
itask = int64(1);
iout  = 1;
xout  = 2;
y     = [1; 0; 0];
itol  = int64(1);
rtol  = [1e-04];
atol  = [1e-07];

% Output header and initial results.
fprintf('\n      x           y(1)       y(2)          y(3)\n');
fprintf(' %8.3f %12.4f %12.2e %11.4f\n', t, y);

% Prepare to store results for plotting.
ncall = 1;
ykeep = y;
tkeep = [t];

% bogus initial entry values for unitialized input parameters.
irevcm = int64(-999); imon = irevcm; inln = irevcm; ires = irevcm;

% Ask for warning messages only.
itrace = int64(0);

% The reverse communication process is controlled by the value of irevcm
% (which is 0 for the final exit).
while (irevcm ~= 0)
  [t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, imon, inln, ires, ...
   irevcm, ifail] = d02nm(...
                          t, tout, y, ydot, rwork, rtol, atol, itol, ...
                          inform, ysave, wkjac, jacpvt, imon, inln, ...
                          ires, irevcm, itask, itrace, 'neq', neq, ...
                          'sdysav', sdysav, 'nwkjac', nwkjac);

  % Evaluate the derivative, then use irevcm to place it.
  f(1) = -0.04*y(1) + 1.0e4*y(2)*y(3);
  f(2) =  0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2);
  f(3) =  3.0e7*y(2)*y(2);
  if (irevcm == 1 || irevcm == 3 )
    lrw = 50+2*neq;
    rwork(lrw+1:lrw+neq) = f;
  elseif (irevcm == 4)
    lrw = 50+neq;
    rwork(lrw+1:lrw+neq) = f;
  elseif (irevcm == 5)
    ydot = f;
  elseif (irevcm == 9 && imon == 1)
    % Extract useful information about the current step.
    tc    = rwork(19);
    hlast = rwork(15);
    hnext = rwork(16);
    nqu   = int64(rwork(10));

    % For low values of the independent variable, store the results
    % for plotting.
    if tc < 2.0
      ncall = ncall + 1;
      tkeep(ncall,1) = tc;
      ykeep(:,ncall) = y;
    end

    % xout is intermediate output point.
    while (xout <= tc)
      lrw = 50+neq;
      [sol, ifail] = d02xk(...
                           xout, neq, ysave, rwork(lrw+1:lrw+neq), ...
                           tc, nqu, hlast, hnext);
      % Output interpolated result, and save it for plotting.
      fprintf(' %8.3f %12.4f %12.2e %11.4f\n', xout, sol);
      ncall = ncall + 1;
      tkeep(ncall) = xout;
      ykeep(:,ncall) = sol;
      iout = iout + 1;
      xout = iout*2;
    end
  end
end

% d02ny is an integrator diagnostic routine which can be called
% after d02nm.  Output results.
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] ...
  = d02ny(...
          neq, neqmax, rwork, inform);
fprintf('\nDiagnostic information\n integration status:\n');
fprintf('   last and next step sizes = %8.5f, %8.5f\n',hu, h);
fprintf('   integration stopped at x = %8.5f\n',tcur);
fprintf(' algorithm statistics:\n');
fprintf('   number of time-steps and Newton iterations  = %5d %5d\n', ...
        nst,niter);
fprintf('   number of residual and jacobian evaluations = %5d %5d\n', ...
        nre,nje);
fprintf('   order of method last used and next to use   = %5d %5d\n', ...
        nqu,nq);
fprintf('   component with largest error                = %5d\n',imxer);
% Plot results.
fig1 = figure;
display_plot(tkeep, ykeep)


function display_plot(x, y)
  % Plot one of the curves and then add the other two.
  hline1 = semilogx(x, y(1,:));
  hold on
  [haxes, hline2, hline3] = plotyy(x,y(3,:),x,y(2,:), ...
                                   'semilogx','semilogx');
  % Set the axis limits and the tick specifications to beautify the plot.
  set(haxes(1), 'YLim', [-0.05 1.3]);
  set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
  set(haxes(1), 'YTick', [0.0:0.2:1.2]);
  set(haxes(2), 'YLim', [0.0 5e-5]);
  set(haxes(2), 'YMinorTick', 'on');
  set(haxes(2), 'YTick', [5e-6:5e-6:3.5e-5]);
  for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 12]);
    set(haxes(iaxis), 'XTick', [0.0001 0.001 0.01 0.1 1 10]);
  end
  set(gca, 'box', 'off');
  % Add title.
  title({'Stiff Robertson Problem (Explicit ODE):', ...
         'BDF with Newton Iterations'},'Position',[0.05,1.1]);
  % Label the x axis, and both y axes.
  xlabel('x');
  ylabel(haxes(1),'Solution (a,c)');
  ylabel(haxes(2),'Solution (b)');
  % Add a legend.
  legend({'a','c','b'},'Position',[0.45,0.45,0.15,0.15]);
  % Set some features of the three lines.
  set(hline1, 'Marker', '+','Linestyle','-','Color','red');
  set(hline2, 'Marker', 'o','Linestyle',':','Color','green');
  set(hline3, 'Marker', 'x','Linestyle','--','Color','blue');
  hold off

d02nm example results


      x           y(1)       y(2)          y(3)
    0.000       1.0000     0.00e+00      0.0000
    2.000       0.9416     2.70e-05      0.0584
    4.000       0.9055     2.24e-05      0.0945
    6.000       0.8793     1.96e-05      0.1207
    8.000       0.8585     1.77e-05      0.1414
   10.000       0.8414     1.62e-05      0.1586

Diagnostic information
 integration status:
   last and next step sizes =  0.90178,  0.90178
   integration stopped at x = 10.76621
 algorithm statistics:
   number of time-steps and Newton iterations  =    55    78
   number of residual and jacobian evaluations =   128    16
   order of method last used and next to use   =     4     4
   component with largest error                =     3