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

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_ivp_stiff_imp_sparjac (d02nj)

## Purpose

nag_ode_ivp_stiff_imp_sparjac (d02nj) is a forward communication function for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations when the Jacobian is a sparse matrix.

## Syntax

[t, tout, y, ydot, rwork, inform, ysav, wkjac, jacpvt, lderiv, ifail] = d02nj(t, tout, y, ydot, rwork, rtol, atol, itol, inform, resid, ysav, jac, wkjac, jacpvt, monitr, lderiv, itask, itrace, 'neq', neq, 'sdysav', sdysav)
[t, tout, y, ydot, rwork, inform, ysav, wkjac, jacpvt, lderiv, ifail] = nag_ode_ivp_stiff_imp_sparjac(t, tout, y, ydot, rwork, rtol, atol, itol, inform, resid, ysav, jac, wkjac, jacpvt, monitr, lderiv, itask, itrace, 'neq', neq, 'sdysav', sdysav)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: nwkjac, njcpvt have been removed from the interface; neq has been made optional
.

## Description

nag_ode_ivp_stiff_imp_sparjac (d02nj) is a general purpose function for integrating the initial value problem for a stiff system of implicit ordinary differential equations coupled with algebraic equations, written in the form
 A(t,y)y′ = g(t,y). $A(t,y)y′=g(t,y).$
It is designed specifically for the case where the resulting Jacobian is a sparse matrix (see the description of jac).
Both interval and step oriented modes of operation are available and also modes designed to permit intermediate output within an interval oriented mode.
An outline of a typical calling program for nag_ode_ivp_stiff_imp_sparjac (d02nj) is given below. It calls the sparse matrix linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu), the Backward Differentiation Formula (BDF) integrator setup function nag_ode_ivp_stiff_bdf (d02nv), its diagnostic counterpart nag_ode_ivp_stiff_integ_diag (d02ny), and the sparse matrix linear algebra diagnostic function nag_ode_ivp_stiff_sparjac_diag (d02nx).
```.
.
.
[...] = d02nv(...);
[...] = d02nu(...);
[..., ifail] = d02nj(...);
if (ifail ~= 1 and ifail < 14)
[...] = d02nx(...);
[...] = d02ny(...);
end
.
.
.
```
The linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu) and one of the integrator setup functions, nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw), must be called prior to the call of nag_ode_ivp_stiff_imp_sparjac (d02nj). Either or both of the integrator diagnostic function nag_ode_ivp_stiff_integ_diag (d02ny), or the sparse matrix linear algebra diagnostic function nag_ode_ivp_stiff_sparjac_diag (d02nx), may be called after the call to nag_ode_ivp_stiff_imp_sparjac (d02nj). There is also a function, nag_ode_ivp_stiff_contin (d02nz), designed to permit you to change step size on a continuation call to nag_ode_ivp_stiff_imp_sparjac (d02nj) without restarting the integration process.

## References

See the D02M–N sub-chapter Introduction.

## Parameters

### Compulsory Input Parameters

1:     t – double scalar
t$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:     tout – double scalar
The next value of t$t$ at which a computed solution is desired. For the initial t$t$, the input value of tout is used to determine the direction of integration. Integration is permitted in either direction (see also itask).
3:     y(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
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:     ydot(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
If lderiv(1) = true${\mathbf{lderiv}}\left(1\right)=\mathbf{true}$, ydot must contain approximations to the time derivatives y${y}^{\prime }$ of the vector y$y$.
If lderiv(1) = false${\mathbf{lderiv}}\left(1\right)=\mathbf{false}$, ydot need not be set on entry.
5:     rwork(50 + 4 × neq$50+4×{\mathbf{neq}}$) – double array
6:     rtol( : $:$) – double array
Note: the dimension of the array rtol must be at least 1$1$ if itol = 1${\mathbf{itol}}=1$ or 2$2$, and at least neq${\mathbf{neq}}$ otherwise.
The relative local error tolerance.
Constraint: rtol(i)0.0${\mathbf{rtol}}\left(i\right)\ge 0.0$ for all relevant i$i$ (see itol).
7:     atol( : $:$) – double array
Note: the dimension of the array atol must be at least 1$1$ if itol = 1${\mathbf{itol}}=1$ or 3$3$, and at least neq${\mathbf{neq}}$ otherwise.
The absolute local error tolerance.
Constraint: atol(i)0.0${\mathbf{atol}}\left(i\right)\ge 0.0$ for all relevant i$i$ (see itol).
8:     itol – int64int32nag_int scalar
A value to indicate the form of the local error test. itol indicates to nag_ode_ivp_stiff_imp_sparjac (d02nj) whether to interpret either or both of rtol or atol as a vector or a scalar. The error test to be satisfied is ei / wi < 1.0$‖{e}_{i}/{w}_{i}‖<1.0$, where wi${w}_{i}$ is defined as follows:
 itol rtol atol wi${w}_{i}$ 1 scalar scalar rtol(1) × |yi| + atol(1)${\mathbf{rtol}}\left(1\right)×|{y}_{i}|+{\mathbf{atol}}\left(1\right)$ 2 scalar vector rtol(1) × |yi| + atol(i)${\mathbf{rtol}}\left(1\right)×|{y}_{i}|+{\mathbf{atol}}\left(i\right)$ 3 vector scalar rtol(i) × |yi| + atol(1)${\mathbf{rtol}}\left(i\right)×|{y}_{i}|+{\mathbf{atol}}\left(1\right)$ 4 vector vector rtol(i) × |yi| + atol(i)${\mathbf{rtol}}\left(i\right)×|{y}_{i}|+{\mathbf{atol}}\left(i\right)$
ei${e}_{i}$ is an estimate of the local error in yi${y}_{i}$, computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup function.
Constraint: itol = 1${\mathbf{itol}}=1$, 2$2$, 3$3$ or 4$4$.
9:     inform(23$23$) – int64int32nag_int array
10:   resid – function handle or string containing name of m-file
resid must evaluate the residual
 r = g(t,y) − A(t,y)y′ $r=g(t,y)-A(t,y)y′$
in one case and
 r = − A(t,y)y′ $r=-A(t,y)y′$
in another.
[r, ires] = resid(neq, t, y, ydot, ires)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of equations being solved.
2:     t – double scalar
t$t$, the current value of the independent variable.
3:     y(neq) – double array
The value of yi${y}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
4:     ydot(neq) – double array
The value of yi${y}_{\mathit{i}}^{\prime }$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, at t$t$.
5:     ires – int64int32nag_int scalar
The form of the residual that must be returned in array r.
ires = -1${\mathbf{ires}}=-1$
The residual defined in equation (2) must be returned.
ires = 1${\mathbf{ires}}=1$
The residual defined in equation (1) must be returned.

Output Parameters

1:     r(neq) – double array
r(i)${\mathbf{r}}\left(\mathit{i}\right)$ must contain the i$\mathit{i}$th component of r$r$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, where
 r = g(t,y) − A(t,y)y′ $r=g(t,y)-A(t,y)y′$ (1)
or
 r = − A(t,y)y′ $r=-A(t,y)y′$ (2)
and where the definition of r$r$ is determined by the input value of ires.
2:     ires – int64int32nag_int scalar
Should be unchanged unless one of the following actions is required of the integrator, in which case ires should be set accordingly.
ires = 2${\mathbf{ires}}=2$
Indicates to the integrator that control should be passed back immediately to the calling (sub)program with the error indicator set to ${\mathbf{ifail}}={\mathbf{11}}$.
ires = 3${\mathbf{ires}}=3$
Indicates to the integrator that an error condition has occurred in the solution vector, its time derivative or in the value of t$t$. The integrator will use a smaller time step to try to avoid this condition. If this is not possible, the integrator returns to the calling (sub)program with the error indicator set to ${\mathbf{ifail}}={\mathbf{7}}$.
ires = 4${\mathbf{ires}}=4$
Indicates to the integrator to stop its current operation and to enter monitr immediately with parameter imon = 2${\mathbf{imon}}=-2$.
11:   ysav(ldysav,sdysav) – double array
12:   jac – function handle or string containing name of m-file
jac must evaluate the Jacobian of the system. If this option is not required, the actual argument for jac must be the string 'd02njz'. (nag_ode_ivp_stiff_imp_sparjac_dummy_jac (d02njz) is included in the NAG Toolbox.) You must indicate to the integrator whether this option is to be used by setting the parameter jceval appropriately in a call to the sparse linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu).
First we must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, y${y}^{\prime }$, generated internally, has the form
 y′ = (y − z) / (hd) , $y′ = (y-z) / (hd) ,$
where h$h$ is the current step size and d$d$ is a parameter that depends on the integration method in use. The vector y$y$ is the current solution and the vector z$z$ depends on information from previous time steps. This means that d/(dy) (​ ​) = (hd) d/(dy) (​ ​) $\frac{d}{d{y}^{\prime }}\left(\text{​ ​}\right)=\left(hd\right)\frac{d}{dy}\left(\text{​ ​}\right)$. The system of nonlinear equations that is solved has the form
 A (t,y) y′ − g (t,y) = 0 $A (t,y) y′ - g (t,y) = 0$
but it is solved in the form
 r (t,y) = 0 , $r (t,y) = 0 ,$
where r$r$ is the function defined by
 r (t,y) = (hd) ( A (t,y) (y − z) / (hd) − g (t,y) ) . $r (t,y) = (hd) ( A (t,y) (y-z) / (hd) - g (t,y) ) .$
It is the Jacobian matrix (r)/(y) $\frac{\partial r}{\partial y}$ that you must supply in jac as follows:
(ri)/(yj) = aij (t,y) + (hd) ()/(yj)
 ( neq ) ∑ aik(t,y)yk ′ − gi(t,y) k = 1
.
$∂ri ∂yj = aij (t,y) + (hd) ∂ ∂yj ( ∑ k=1 neq aik (t,y) yk′ - gi (t,y) ) .$
[pdj] = jac(neq, t, y, ydot, h, d, j, pdj)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of equations being solved.
2:     t – double scalar
t$t$, the current value of the independent variable.
3:     y(neq) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4:     ydot(neq) – double array
The derivative of the solution at the current point t$t$.
5:     h – double scalar
The current step size.
6:     d – double scalar
The parameter d$d$ which depends on the integration method.
7:     j – int64int32nag_int scalar
The column of the Jacobian that jac must return in the array pdj.
8:     pdj(neq) – double array
Is set to zero.

Output Parameters

1:     pdj(neq) – double array
pdj(i)${\mathbf{pdj}}\left(i\right)$ should be set to the (i,j)$\left(i,j\right)$th element of the Jacobian, where j$j$ is given by j. Only nonzero elements of this array need be set, since it is preset to zero before the call to jac.
13:   wkjac(nwkjac) – double array
The actual size depends on whether the sparsity structure is supplied or whether it is to be estimated. An appropriate value for nwkjac is described in the specification of the linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu). This value must be the same as that supplied to nag_ode_ivp_stiff_sparjac_setup (d02nu).
14:   jacpvt(njcpvt) – int64int32nag_int array
The actual size depends on whether the sparsity structure is supplied or whether it is to be estimated. An appropriate value for njcpvt is described in the specification for the linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu). This value must be same as that supplied to nag_ode_ivp_stiff_sparjac_setup (d02nu).
15:   monitr – function handle or string containing name of m-file
monitr performs tasks requested by you. If this option is not required, then the actual argument for monitr must be the string 'd02nby'. (nag_ode_ivp_stiff_exp_fulljac_dummy_monit (d02nby) is included in the NAG Toolbox.)
[hnext, y, imon, inln, hmin, hmax] = monitr(neq, ldysav, t, hlast, hnext, y, ydot, ysav, r, acor, imon, hmin, hmax, nqu)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of equations being solved.
2:     ldysav – int64int32nag_int scalar
An upper bound on the number of equations to be solved.
3:     t – double scalar
The current value of the independent variable.
4:     hlast – double scalar
The last step size successfully used by the integrator.
5:     hnext – double scalar
The step size that the integrator proposes to take on the next step.
6:     y(neq) – double array
y$y$, the values of the dependent variables evaluated at t$t$.
7:     ydot(neq) – double array
The time derivatives y${y}^{\prime }$ of the vector y$y$.
8:     ysav(ldysav,sdysav$\mathit{sdysav}$) – double array
Workspace to enable you to carry out interpolation using either of the functions nag_ode_ivp_stiff_nat_interp (d02xj) or nag_ode_ivp_stiff_c1_interp (d02xk).
9:     r(neq) – double array
If imon = 0${\mathbf{imon}}=0$ and inln = 3${\mathbf{inln}}=3$, then the first neq elements contain the residual vector A(t,y)yg(t,y)$A\left(t,y\right){y}^{\prime }-g\left(t,y\right)$.
10:   acor(neq,2$2$) – double array
With imon = 1${\mathbf{imon}}=1$, acor(i,1)${\mathbf{acor}}\left(i,1\right)$ contains the weight used for the i$i$th equation when the norm is evaluated, and acor(i,2)${\mathbf{acor}}\left(i,2\right)$ contains the estimated local error for the i$i$th equation. The scaled local error at the end of a timestep may be obtained by calling the double function nag_ode_ivp_stiff_errest (d02za) as follows:
` [errloc, ifail] = d02za(acor(1:neq,2), acor(1:neq,1)); % Check ifail before proceeding `
11:   imon – int64int32nag_int scalar
A flag indicating under what circumstances monitr was called:
imon = -2${\mathbf{imon}}=-2$
Entry from the integrator after ires = 4${\mathbf{ires}}=4$ (set in resid) caused an early termination (this facility could be used to locate discontinuities).
imon = -1${\mathbf{imon}}=-1$
The current step failed repeatedly.
imon = 0${\mathbf{imon}}=0$
Entry after a call to the internal nonlinear equation solver (see inln).
imon = 1${\mathbf{imon}}=1$
The current step was successful.
12:   hmin – double scalar
The minimum step size to be taken on the next step.
13:   hmax – double scalar
The maximum step size to be taken on the next step.
14:   nqu – int64int32nag_int scalar
The order of the integrator used on the last step. This is supplied to enable you to carry out interpolation using either of the functions nag_ode_ivp_stiff_nat_interp (d02xj) or nag_ode_ivp_stiff_c1_interp (d02xk).

Output Parameters

1:     hnext – double scalar
The next step size to be used. If this is different from the input value, then imon must be set to 4$4$.
2:     y(neq) – double array
These values must not be changed unless imon is set to 2$2$.
3:     imon – int64int32nag_int scalar
May be reset to determine subsequent action in nag_ode_ivp_stiff_imp_sparjac (d02nj).
imon = -2${\mathbf{imon}}=-2$
Integration is to be halted. A return will be made from the integrator to the calling (sub)program with ${\mathbf{ifail}}={\mathbf{12}}$.
imon = -1${\mathbf{imon}}=-1$
Allow the integrator to continue with its own internal strategy. The integrator will try up to three restarts unless imon-1${\mathbf{imon}}\ne -1$ on exit.
imon = 0${\mathbf{imon}}=0$
Return to the internal nonlinear equation solver, where the action taken is determined by the value of inln (see inln).
imon = 1${\mathbf{imon}}=1$
Normal exit to the integrator to continue integration.
imon = 2${\mathbf{imon}}=2$
Restart the integration at the current time point. The integrator will restart from order 1$1$ when this option is used. The solution y, provided by monitr, will be used for the initial conditions.
imon = 3${\mathbf{imon}}=3$
Try to continue with the same step size and order as was to be used before the call to monitr. hmin and hmax may be altered if desired.
imon = 4${\mathbf{imon}}=4$
Continue the integration but using a new value of hnext and possibly new values of hmin and hmax.
4:     inln – int64int32nag_int scalar
The action to be taken by the internal nonlinear equation solver when monitr is exited with imon = 0${\mathbf{imon}}=0$. By setting inln = 3${\mathbf{inln}}=3$ and returning to the integrator, the residual vector is evaluated and placed in the array r, and then monitr is called again. At present this is the only option available: inln must not be set to any other value.
5:     hmin – double scalar
The minimum step size to be used. If this is different from the input value, then imon must be set to 3$3$ or 4$4$.
6:     hmax – double scalar
The maximum step size to be used. If this is different from the input value, then imon must be set to 3$3$ or 4$4$. If hmax is set to zero, no limit is assumed.
16:   lderiv(2$2$) – logical array
lderiv(1)${\mathbf{lderiv}}\left(1\right)$ must be set to true if you have supplied both an initial y$y$ and an initial y${y}^{\prime }$. lderiv(1)${\mathbf{lderiv}}\left(1\right)$ must be set to false if only the initial y$y$ has been supplied.
lderiv(2)${\mathbf{lderiv}}\left(2\right)$ must be set to true if the integrator is to use a modified Newton method to evaluate the initial y$y$ and y${y}^{\prime }$. Note that y$y$ and y${y}^{\prime }$, if supplied, are used as initial estimates. This method involves taking a small step at the start of the integration, and if itask = 6${\mathbf{itask}}=6$ on entry, t and tout will be set to the result of taking this small step. lderiv(2)${\mathbf{lderiv}}\left(2\right)$ must be set to false if the integrator is to use functional iteration to evaluate the initial y$y$ and y${y}^{\prime }$, and if this fails a modified Newton method will then be attempted. lderiv(2) = true${\mathbf{lderiv}}\left(2\right)=\mathbf{true}$ is recommended if there are implicit equations or the initial y$y$ and y${y}^{\prime }$ are zero.
17:   itask – int64int32nag_int scalar
The task to be performed by the integrator.
itask = 1${\mathbf{itask}}=1$
Normal computation of output values of y(t)$y\left(t\right)$ at t = tout$t={\mathbf{tout}}$ (by overshooting and interpolating).
itask = 2${\mathbf{itask}}=2$
Take one step only and return.
itask = 3${\mathbf{itask}}=3$
Stop at the first internal integration point at or beyond t = tout$t={\mathbf{tout}}$ and return.
itask = 4${\mathbf{itask}}=4$
Normal computation of output values of y(t)$y\left(t\right)$ at t = tout$t={\mathbf{tout}}$ but without overshooting t = tcrit$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 may be equal to or beyond tout, but not before it, in the direction of integration.
itask = 5${\mathbf{itask}}=5$
Take one step only and return, without passing tcrit. tcrit must be specified as under itask = 4${\mathbf{itask}}=4$.
itask = 6${\mathbf{itask}}=6$
The integrator will solve for the initial values of y$y$ and y${y}^{\prime }$ only and then return to the calling (sub)program without doing the integration. This option can be used to check the initial values of y$y$ and y${y}^{\prime }$. Functional iteration or a ‘small’ backward Euler method used in conjunction with a damped Newton iteration is used to calculate these values (see lderiv). Note that if a backward Euler step is used then the value of t$t$ will have been advanced a short distance from the initial point.
Note:  if nag_ode_ivp_stiff_imp_sparjac (d02nj) is recalled with a different value of itask (and tout altered), then the initialization procedure is repeated, possibly leading to different initial conditions.
Constraint: 1itask6$1\le {\mathbf{itask}}\le 6$.
18:   itrace – int64int32nag_int scalar
The level of output that is printed by the integrator. itrace may take the value 1$-1$, 0$0$, 1$1$, 2$2$ or 3$3$.
itrace < 1${\mathbf{itrace}}<-1$
1$-1$ is assumed and similarly if itrace > 3${\mathbf{itrace}}>3$, then 3$3$ is assumed.
itrace = 1${\mathbf{itrace}}=-1$
No output is generated.
itrace = 0${\mathbf{itrace}}=0$
Only warning messages are printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
itrace > 0${\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:     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.)
The number of differential equations to be solved.
Constraint: neq1${\mathbf{neq}}\ge 1$.
2:     sdysav – int64int32nag_int scalar
Default: The second dimension of the array ysav.
The second dimension of the array ysav as declared in the (sub)program from which nag_ode_ivp_stiff_imp_sparjac (d02nj) is called. An appropriate value for sdysav is described in the specifications of the integrator setup functions nag_ode_ivp_stiff_dassl (d02mv), 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.

### Input Parameters Omitted from the MATLAB Interface

ldysav nwkjac njcpvt

### Output Parameters

1:     t – double scalar
The value at which the computed solution y$y$ is returned (usually at tout).
2:     tout – double scalar
Normally unchanged. However, when itask = 6${\mathbf{itask}}=6$, then tout contains the value of t at which initial values have been computed without performing any integration. See descriptions of itask and lderiv.
3:     y(neq) – double array
The computed solution vector, evaluated at t (usually ${\mathbf{t}}={\mathbf{tout}}$).
4:     ydot(neq) – double array
The time derivatives y${y}^{\prime }$ of the vector y$y$ at the last integration point.
5:     rwork(50 + 4 × neq$50+4×{\mathbf{neq}}$) – double array
6:     inform(23$23$) – int64int32nag_int array
7:     ysav(ldysav,sdysav) – double array
8:     wkjac(nwkjac) – double array
Communication array, used to store information between calls to nag_ode_ivp_stiff_imp_sparjac (d02nj).
9:     jacpvt(njcpvt) – int64int32nag_int array
Communication array, used to store information between calls to nag_ode_ivp_stiff_imp_sparjac (d02nj).
10:   lderiv(2$2$) – logical array
lderiv(1)${\mathbf{lderiv}}\left(1\right)$ is normally unchanged. However if itask = 6${\mathbf{itask}}=6$ and internal initialization was successful then lderiv(1) = true${\mathbf{lderiv}}\left(1\right)=\mathbf{true}$.
lderiv(2) = true${\mathbf{lderiv}}\left(2\right)=\mathbf{true}$, if implicit equations were detected. Otherwise lderiv(2) = false${\mathbf{lderiv}}\left(2\right)=\mathbf{false}$.
11:   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:

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

ifail = 1${\mathbf{ifail}}=1$
An illegal input was detected on entry, or after an internal call to monitr. If itrace > 1${\mathbf{itrace}}>-1$, then the form of the error will be detailed on the current error message unit (see nag_file_set_unit_error (x04aa)).
ifail = 2${\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 ifail = 3${\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 y(1),y(2),,y(neq)${\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 ifail = 4${\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 ifail = 5${\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 ifail = 6${\mathbf{ifail}}=6$
Some error weight wi${w}_{i}$ became zero during the integration (see the description of itol). Pure relative error control (atol(i) = 0.0${\mathbf{atol}}\left(i\right)=0.0$) was requested on a variable (the i$i$th) which has now vanished. The integration was successful as far as t.
ifail = 7${\mathbf{ifail}}=7$
resid set its error flag (ires = 3${\mathbf{ires}}=3$) continually despite repeated attempts by the integrator to avoid this.
ifail = 8${\mathbf{ifail}}=8$
lderiv(1) = false${\mathbf{lderiv}}\left(1\right)=\mathbf{false}$ on entry but the internal initialization function was unable to initialize y${y}^{\prime }$ (more detailed information may be directed to the current error message unit, see nag_file_set_unit_error (x04aa)).
ifail = 9${\mathbf{ifail}}=9$
A singular Jacobian (r)/(y) $\frac{\partial r}{\partial y}$ has been encountered. You should check the problem formulation and Jacobian calculation.
ifail = 10${\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 ifail = 11${\mathbf{ifail}}=11$
resid signalled the integrator to halt the integration and return (ires = 2${\mathbf{ires}}=2$). Integration was successful as far as t.
W ifail = 12${\mathbf{ifail}}=12$
monitr set imon = 2${\mathbf{imon}}=-2$ and so forced a return but the integration was successful as far as t.
W ifail = 13${\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 itask2${\mathbf{itask}}\ne 2$ or 5$5$.)
ifail = 14${\mathbf{ifail}}=14$
The values of rtol and atol are so small that nag_ode_ivp_stiff_imp_sparjac (d02nj) is unable to start the integration.
ifail = 15${\mathbf{ifail}}=15$
The linear algebra setup function nag_ode_ivp_stiff_sparjac_setup (d02nu) was not called before the call to nag_ode_ivp_stiff_imp_sparjac (d02nj).

## Accuracy

The accuracy of the numerical solution may be controlled by a careful choice of the parameters rtol and atol, and to a much lesser extent by the choice of norm. You are advised to use scalar error control unless the components of the solution are expected to be poorly scaled. For the type of decaying solution typical of many stiff problems, relative error control with a small absolute error threshold will be most appropriate (that is, you are advised to choose itol = 1${\mathbf{itol}}=1$ with atol(1)${\mathbf{atol}}\left(1\right)$ small but positive).

Since numerical stability and memory are often conflicting requirements when solving ordinary differential systems where the Jacobian matrix is sparse we provide a diagnostic function, nag_ode_ivp_stiff_sparjac_diag (d02nx), whose aim is to inform you how much memory is required to solve the problem and to give you some indicators of numerical stability.
In general, you are advised to choose the BDF option (setup function nag_ode_ivp_stiff_bdf (d02nv)) but if efficiency is of great importance and especially if it is suspected that ()/(y) (A1g) $\frac{\partial }{\partial y}\left({A}^{-1}g\right)$ has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup function nag_ode_ivp_stiff_blend (d02nw)).

## Example

This example solves the well-known stiff Robertson problem written as a mixed differential/algebraic system in implicit form
 r1 = a + b + c − 1.0 r2 = 0.04a − 1.0E4bc − 3.0E7b2 − b′ r3 = 0.04a − 1.0E4bc − 3.0E7b2 − c′
$r1 = a+b+c-1.0 r2 = 0.04a-1.0E4bc-3.0E7⁢b2-b′ r3 = 0.04a-1.0E4bc-3.0E7⁢b2-c′$
exploiting the fact that, from the initial conditions a = 1.0$a=1.0$ and b = c = 0.0$b=c=0.0$, we know that a + b + c = 1$a+b+c=1$ for all time. We integrate over the range [0,10.0]$\left[0,10.0\right]$ with vector relative error control and scalar absolute error control (itol = 3${\mathbf{itol}}=3$) and using the BDF integrator (setup function nag_ode_ivp_stiff_bdf (d02nv)) and a modified Newton method. The Jacobian is evaluated, in turn, using the 'A' (Analytical) and 'F' (Full information) options. We provide a monitor function to terminate the integration when the value of the component a falls below 0.9$0.9$.
```function nag_ode_ivp_stiff_imp_sparjac_example
% For communication with monitr.
global ncall tkeep ykeep;

neq = 3;
neqmax = neq;
nrw = 50+4*neqmax;
ninf = 23;
nelts = 8;
njcpvt = 150;
nwkjac = 100;
nia = neq+1;
nja = nelts;
maxord = 5;
sdysav = maxord+1;
maxstp = 200;
mxhnil = 5;
ldysav = neq;

h0 = 0.0;
hmax = 10.0;
hmin = 1.0e-10;
tcrit = 0.0;
petzld = true;
eta = 1.0e-4;
u = 0.1;
sens = 1.0e-6;
lblock = true;

% Initialize variables and arrays.
const = zeros(6, 1);
rwork = zeros(nrw, 1);
wkjac = zeros(nwkjac, 1);
ydot = zeros(neq, 1);
ysave = zeros(ldysav, sdysav);
inform = zeros(ninf, 1);
jacpvt = zeros(njcpvt, 1);
algequ = zeros(neq, 1);
ia = [1; 3; 6; 9];
ja = [1; 2; 1; 2; 3; 1; 2; 3];

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

% We integrate to tout by overshooting (itask=1) using the backward
% differentiation formula (BDF) formulae with a Newton method.  Also, we
% set petzld to true so that the Petzold error test is used (because an
% algebraic equation is defined in the system).  Default values for the
% array const are used.  Employ vector relative tolerance and scalar
% absolute tolerance.  The monitr routine is used to force a return
% when the first component of the system falls below the value 0.9.

% The Jacobian is supplied by jac.  In the first case, its sparsity
% structure is determined internally by calls to that function.  In the
% second, it's given explictly (by the arrays ia and ja).  This is
% specified by the value of jceval in the call to nag_ode_ivp_stiff_sparjac_setup.
jceval = [{'Analytical'}, {'Full information'}];
title = [{'Analytic Jacobian, structure not supplied\n\n'},...
{'Analytic Jacobian, structure supplied\n\n'}];

for icase = 1:2
% Set initial parameter values.
t = 0;
tout = 10.0;
y = [1; 0; 0];
lderiv = [false; false];
itol = 3;
rtol = [0.0001; 0.001; 0.0001];
atol = [1e-07];
isplit = 0;

% nag_ode_ivp_stiff_bdf is a setup routine (specifically for the backward
% differentiation formula (BDF) integrator) to be called prior to nag_ode_ivp_stiff_imp_sparjac.
[const, rwork, ifail] = nag_ode_ivp_stiff_bdf(int64(neqmax), int64(sdysav), ...
int64(maxord), 'Newton', petzld, const, tcrit, hmin, ...
hmax, h0, int64(maxstp), int64(mxhnil), 'Average-L2', rwork);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: nag_ode_ivp_stiff_bdf returned with ifail = %1d ',ifail);
end

% nag_ode_ivp_stiff_sparjac_setup is a setup routine (specifically for sparse matrix linear
% algebra) to be called prior to nag_ode_ivp_stiff_imp_sparjac.
[jacpvt, rwork, ifail] = nag_ode_ivp_stiff_sparjac_setup(int64(neq), int64(neqmax), ...
jceval{icase}, int64(nwkjac), int64(ia), int64(ja), ...
int64(njcpvt), sens, u, eta, lblock, rwork);

% Output header and initial results.
fprintf(title{icase});
fprintf('    X           Y(1)           Y(2)           Y(3)\n');
fprintf('%8.3f       ',t);
fprintf('%1.5f       ',y);
fprintf('\n');

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

% Set itrace to ask for only warning messages from the integrator.
% Also, turn off warnings from NAG routines, because this call is set
% up to be terminated by the monitor routine once integration has
% reached a=0.9 (this means nag_ode_ivp_stiff_imp_sparjac will terminate with ifail=12).
itrace = 0;
warning('off', 'NAG:warning');
[tOut, toutOut, yOut, ydotOut, rworkOut, informOut, ysaveOut, ...
wkjacOut, jacpvtOut, lderivOut, ifail] = nag_ode_ivp_stiff_imp_sparjac(t, tout, y, ...
ydot, rwork, rtol, atol, int64(itol), int64(inform), ...
@resid, ysave, @jac, wkjac, int64(jacpvt), ...
@monitr, lderiv, int64(itask), int64(itrace));
warning('on', 'NAG:warning');

if (ifail == 0 || ifail == 12)
% Output results.
fprintf('%8.3f       ',tOut);
fprintf('%1.5f       ',yOut);
fprintf('\n');

% nag_ode_ivp_stiff_integ_diag is an integrator diagnostic routine which can be called
% after nag_ode_ivp_stiff_imp_sparjac.
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, ...
algequ, ifail] = nag_ode_ivp_stiff_integ_diag(int64(neq), int64(neqmax), ...
rworkOut, informOut);

fprintf('\n hused = %1.5e   hnext = %1.5e   tcur = %1.5e\n',...
hu,h,tcur);
fprintf(' nst = %1.0f   nre = %1.0f   nje = %1.0f\n',nst,nre,nje);
fprintf(' nqu = %1.0f   nq = %1.0f   niter = %1.0f\n',nqu,nq,niter);
fprintf(' Max Err Comp = %1.0f  \n\n',imxer);

% nag_ode_ivp_stiff_sparjac_diag is a sparse linear algebra diagnostic routine which can be
% called after nag_ode_ivp_stiff_imp_sparjac.
icall = 0;
[liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow,...
nblock] = nag_ode_ivp_stiff_sparjac_diag(int64(icall), lblock, informOut);

fprintf(' njcpvt (required %1.0f used %1.0f)\n',liwreq,liwusd);
fprintf(' nwkjac (required %1.0f used %1.0f)\n',lrwreq,lrwusd);
fprintf(' No. of LU-decomps %1.0f No. of nonzeros %1.0f\n',nlu,nz);
fprintf([' No. of FCN calls to form Jacobian %1.0f ', ...
'Try isplit %1.0f\n'], ngp,isplit);
fprintf(' Growth est %1.0f No. of blocks on diagonal %1.0f\n\n',...
igrow,nblock);

elseif (ifail == 10 )
icall = 1;
[liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow,...
nblock] = nag_ode_ivp_stiff_sparjac_diag(int64(icall), lblock, informOut);
fprintf(' njcpvt (required %1.0f used %1.0f)\n',liwreq,liwusd);
fprintf(' nwkjac (required %1.0f used %1.0f)\n',lrwreq,lrwusd);

elseif (ifail ~= 0)
fprintf(['Warning: nag_ode_ivp_stiff_imp_sparjac returned with ifail = %1d and ', ...
'T = %1.0f\n\n'], ifail, tOut);
end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function pdj = jac(neq, t, y, ydot, h, d, j, pdj)
% Evaluate the Jacobian.
pdj = zeros(neq,1);
hxd = h*d;
if j == 1
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 0.0 - hxd*(0.04);
elseif j == 2
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 1.0 - hxd*(-1.0e4*y(3)-6.0e7*y(2));
pdj(3) = 0.0 - hxd*(6.0e7*y(2));
elseif j == 3
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 0.0 - hxd*(-1.0e4*y(2));
pdj(3) = 1.0 - hxd*(0.0);
end
function [hnextOut, yOut, imonOut, inln, hminOut, hmaxOut] = ...
monitr(neq, neqmax, t, hlast, hnext, y, ydot, ysave, ...
r, acor, imon, hmin, hmax, nqu)
% For communication with main routine.
global ncall tkeep ykeep;

% If the first component of the solution is less than 0.9, then
% indicate that the integration is to be terminated.
if (y(1) <= 0.9)
imonOut = int64(-2);
else
imonOut = imon;
end

% We have to assign values to all the output arguments (whether
% they're used or not).
inln = int64(0);
hnextOut = hnext;
yOut = y;
hminOut = hmin;
hmaxOut = hmax;

% Store current results for plotting later on.
tkeep(ncall,1) = t;
ykeep(ncall,:) = y;
ncall = ncall + 1;
function [r, ires] = resid(neq, t, y, ydot, ires)
r = zeros(neq,1);
r(1) = 0.0;
r(2) = -ydot(2);
r(3) = -ydot(3);
if ires == 1
r(1) = y(1) + y(2) + y(3) - 1.0d0 + r(1);
r(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) + r(2);
r(3) = 3.0e7*y(2)*y(2) + r(3);
end
function display_plot(x, y)
% 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.
% 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.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.2 0.4 0.6 0.8 1.0 1.2]);
set(haxes(2), 'YLim', [0.0 4e-5]); % xrange must be the same as haxes(1).
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5.0e-6 1.0e-5 1.5e-5 2.0e-5 2.5e-5 3.0e-5 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]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
title(['Stiff Robertson Problem (Implicit DAE): BDF with Newton ', ...
'Iterations'], titleFmt{:});
% Label the x axis, and both y axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution (a,c)', labFmt{:});
ylabel(haxes(2),'Solution (b)', labFmt{:});
% Add a legend.
legend('a','c','b','Location','Best');
% Set some features of the three lines.
set(hline1, 'Marker', '+','Linestyle','-');
set(hline2, 'Marker', '*','Linestyle',':');
set(hline3, 'Marker', 'x','Linestyle','--');
```
```
nag_ode_ivp_stiff_imp_sparjac example program results

Analytic Jacobian, structure not supplied

X           Y(1)           Y(2)           Y(3)
0.000       1.00000       0.00000       0.00000
4.957       0.89208       0.00002       0.10790

hused = 5.99710e-01   hnext = 5.99710e-01   tcur = 4.95659e+00
nst = 52   nre = 132   nje = 12
nqu = 4   nq = 4   niter = 117
Max Err Comp = 3

njcpvt (required 99 used 150)
nwkjac (required 31 used 75)
No. of LU-decomps 12 No. of nonzeros 8
No. of FCN calls to form Jacobian 0 Try isplit 73
Growth est 1034 No. of blocks on diagonal 1

Analytic Jacobian, structure supplied

X           Y(1)           Y(2)           Y(3)
0.000       1.00000       0.00000       0.00000
4.957       0.89208       0.00002       0.10790

hused = 5.99710e-01   hnext = 5.99710e-01   tcur = 4.95659e+00
nst = 52   nre = 131   nje = 12
nqu = 4   nq = 4   niter = 117
Max Err Comp = 3

njcpvt (required 99 used 150)
nwkjac (required 31 used 75)
No. of LU-decomps 12 No. of nonzeros 8
No. of FCN calls to form Jacobian 0 Try isplit 73
Growth est 1034 No. of blocks on diagonal 1

```
```function d02nj_example
% For communication with monitr.
global ncall tkeep ykeep;

neq = 3;
neqmax = neq;
nrw = 50+4*neqmax;
ninf = 23;
nelts = 8;
njcpvt = 150;
nwkjac = 100;
nia = neq+1;
nja = nelts;
maxord = 5;
sdysav = maxord+1;
maxstp = 200;
mxhnil = 5;
ldysav = neq;

h0 = 0.0;
hmax = 10.0;
hmin = 1.0e-10;
tcrit = 0.0;
petzld = true;
eta = 1.0e-4;
u = 0.1;
sens = 1.0e-6;
lblock = true;

% Initialize variables and arrays.
const = zeros(6, 1);
rwork = zeros(nrw, 1);
wkjac = zeros(nwkjac, 1);
ydot = zeros(neq, 1);
ysave = zeros(ldysav, sdysav);
inform = zeros(ninf, 1);
jacpvt = zeros(njcpvt, 1);
algequ = zeros(neq, 1);
ia = [1; 3; 6; 9];
ja = [1; 2; 1; 2; 3; 1; 2; 3];

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

% We integrate to tout by overshooting (itask=1) using the backward
% differentiation formula (BDF) formulae with a Newton method.  Also, we
% set petzld to true so that the Petzold error test is used (because an
% algebraic equation is defined in the system).  Default values for the
% array const are used.  Employ vector relative tolerance and scalar
% absolute tolerance.  The monitr routine is used to force a return
% when the first component of the system falls below the value 0.9.

% The Jacobian is supplied by jac.  In the first case, its sparsity
% structure is determined internally by calls to that function.  In the
% second, it's given explictly (by the arrays ia and ja).  This is
% specified by the value of jceval in the call to d02nu.
jceval = [{'Analytical'}, {'Full information'}];
title = [{'Analytic Jacobian, structure not supplied\n\n'},...
{'Analytic Jacobian, structure supplied\n\n'}];

for icase = 1:2
% Set initial parameter values.
t = 0;
tout = 10.0;
y = [1; 0; 0];
lderiv = [false; false];
itol = 3;
rtol = [0.0001; 0.001; 0.0001];
atol = [1e-07];
isplit = 0;

% d02nv is a setup routine (specifically for the backward
% differentiation formula (BDF) integrator) to be called prior to d02nj.
[const, rwork, ifail] = d02nv(int64(neqmax), int64(sdysav), ...
int64(maxord), 'Newton', petzld, const, tcrit, hmin, ...
hmax, h0, int64(maxstp), int64(mxhnil), 'Average-L2', rwork);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: d02nv returned with ifail = %1d ',ifail);
end

% d02nu is a setup routine (specifically for sparse matrix linear
% algebra) to be called prior to d02nj.
[jacpvt, rwork, ifail] = d02nu(int64(neq), int64(neqmax), ...
jceval{icase}, int64(nwkjac), int64(ia), int64(ja), ...
int64(njcpvt), sens, u, eta, lblock, rwork);

% Output header and initial results.
fprintf(title{icase});
fprintf('    X           Y(1)           Y(2)           Y(3)\n');
fprintf('%8.3f       ',t);
fprintf('%1.5f       ',y);
fprintf('\n');

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

% Set itrace to ask for only warning messages from the integrator.
% Also, turn off warnings from NAG routines, because this call is set
% up to be terminated by the monitor routine once integration has
% reached a=0.9 (this means d02nj will terminate with ifail=12).
itrace = 0;
warning('off', 'NAG:warning');
[tOut, toutOut, yOut, ydotOut, rworkOut, informOut, ysaveOut, ...
wkjacOut, jacpvtOut, lderivOut, ifail] = d02nj(t, tout, y, ...
ydot, rwork, rtol, atol, int64(itol), int64(inform), ...
@resid, ysave, @jac, wkjac, int64(jacpvt), ...
@monitr, lderiv, int64(itask), int64(itrace));
warning('on', 'NAG:warning');

if (ifail == 0 || ifail == 12)
% Output results.
fprintf('%8.3f       ',tOut);
fprintf('%1.5f       ',yOut);
fprintf('\n');

% d02ny is an integrator diagnostic routine which can be called
% after d02nj.
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, ...
algequ, ifail] = d02ny(int64(neq), int64(neqmax), ...
rworkOut, informOut);

fprintf('\n hused = %1.5e   hnext = %1.5e   tcur = %1.5e\n',...
hu,h,tcur);
fprintf(' nst = %1.0f   nre = %1.0f   nje = %1.0f\n',nst,nre,nje);
fprintf(' nqu = %1.0f   nq = %1.0f   niter = %1.0f\n',nqu,nq,niter);
fprintf(' Max Err Comp = %1.0f  \n\n',imxer);

% d02nx is a sparse linear algebra diagnostic routine which can be
% called after d02nj.
icall = 0;
[liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow,...
nblock] = d02nx(int64(icall), lblock, informOut);

fprintf(' njcpvt (required %1.0f used %1.0f)\n',liwreq,liwusd);
fprintf(' nwkjac (required %1.0f used %1.0f)\n',lrwreq,lrwusd);
fprintf(' No. of LU-decomps %1.0f No. of nonzeros %1.0f\n',nlu,nz);
fprintf([' No. of FCN calls to form Jacobian %1.0f ', ...
'Try isplit %1.0f\n'], ngp,isplit);
fprintf(' Growth est %1.0f No. of blocks on diagonal %1.0f\n\n',...
igrow,nblock);

elseif (ifail == 10 )
icall = 1;
[liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow,...
nblock] = d02nx(int64(icall), lblock, informOut);
fprintf(' njcpvt (required %1.0f used %1.0f)\n',liwreq,liwusd);
fprintf(' nwkjac (required %1.0f used %1.0f)\n',lrwreq,lrwusd);

elseif (ifail ~= 0)
fprintf(['Warning: d02nj returned with ifail = %1d and ', ...
'T = %1.0f\n\n'], ifail, tOut);
end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function pdj = jac(neq, t, y, ydot, h, d, j, pdj)
% Evaluate the Jacobian.
pdj = zeros(neq,1);
hxd = h*d;
if j == 1
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 0.0 - hxd*(0.04);
elseif j == 2
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 1.0 - hxd*(-1.0e4*y(3)-6.0e7*y(2));
pdj(3) = 0.0 - hxd*(6.0e7*y(2));
elseif j == 3
pdj(1) = 0.0 - hxd*(1.0);
pdj(2) = 0.0 - hxd*(-1.0e4*y(2));
pdj(3) = 1.0 - hxd*(0.0);
end
function [hnextOut, yOut, imonOut, inln, hminOut, hmaxOut] = ...
monitr(neq, neqmax, t, hlast, hnext, y, ydot, ysave, ...
r, acor, imon, hmin, hmax, nqu)
% For communication with main routine.
global ncall tkeep ykeep;

% If the first component of the solution is less than 0.9, then
% indicate that the integration is to be terminated.
if (y(1) <= 0.9)
imonOut = int64(-2);
else
imonOut = imon;
end

% We have to assign values to all the output arguments (whether
% they're used or not).
inln = int64(0);
hnextOut = hnext;
yOut = y;
hminOut = hmin;
hmaxOut = hmax;

% Store current results for plotting later on.
tkeep(ncall,1) = t;
ykeep(ncall,:) = y;
ncall = ncall + 1;
function [r, ires] = resid(neq, t, y, ydot, ires)
r = zeros(neq,1);
r(1) = 0.0;
r(2) = -ydot(2);
r(3) = -ydot(3);
if ires == 1
r(1) = y(1) + y(2) + y(3) - 1.0d0 + r(1);
r(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) + r(2);
r(3) = 3.0e7*y(2)*y(2) + r(3);
end
function display_plot(x, y)
% 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.
% 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.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.2 0.4 0.6 0.8 1.0 1.2]);
set(haxes(2), 'YLim', [0.0 4e-5]); % xrange must be the same as haxes(1).
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5.0e-6 1.0e-5 1.5e-5 2.0e-5 2.5e-5 3.0e-5 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]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
title(['Stiff Robertson Problem (Implicit DAE): BDF with Newton ', ...
'Iterations'], titleFmt{:});
% Label the x axis, and both y axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution (a,c)', labFmt{:});
ylabel(haxes(2),'Solution (b)', labFmt{:});
% Add a legend.
legend('a','c','b','Location','Best');
% Set some features of the three lines.
set(hline1, 'Marker', '+','Linestyle','-');
set(hline2, 'Marker', '*','Linestyle',':');
set(hline3, 'Marker', 'x','Linestyle','--');
```
```
d02nj example program results

Analytic Jacobian, structure not supplied

X           Y(1)           Y(2)           Y(3)
0.000       1.00000       0.00000       0.00000
4.957       0.89208       0.00002       0.10790

hused = 5.99710e-01   hnext = 5.99710e-01   tcur = 4.95659e+00
nst = 52   nre = 132   nje = 12
nqu = 4   nq = 4   niter = 117
Max Err Comp = 3

njcpvt (required 99 used 150)
nwkjac (required 31 used 75)
No. of LU-decomps 12 No. of nonzeros 8
No. of FCN calls to form Jacobian 0 Try isplit 73
Growth est 1034 No. of blocks on diagonal 1

Analytic Jacobian, structure supplied

X           Y(1)           Y(2)           Y(3)
0.000       1.00000       0.00000       0.00000
4.957       0.89208       0.00002       0.10790

hused = 5.99710e-01   hnext = 5.99710e-01   tcur = 4.95659e+00
nst = 52   nre = 131   nje = 12
nqu = 4   nq = 4   niter = 117
Max Err Comp = 3

njcpvt (required 99 used 150)
nwkjac (required 31 used 75)
No. of LU-decomps 12 No. of nonzeros 8
No. of FCN calls to form Jacobian 0 Try isplit 73
Growth est 1034 No. of blocks on diagonal 1

```

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013