NAG Library Routine Document

D02NHF

1  Purpose

2  Specification

3  Description

!     Declarations

      EXTERNAL RESID, JAC, MONITR
          .
          .
          .
      IFAIL = 0
      CALL D02NVF(...,IFAIL)
      CALL D02NTF(NEQ, NEQMAX, JCEVAL, ML, MU, NWKJAC, NJCPVT, &
      	          RWORK, IFAIL)
      IFAIL = -1
      CALL D02NHF(NEQ, NEQMAX, T, TOUT, Y, YDOT, RWORK, RTOL,  &
                  ATOL, ITOL, INFORM, RESID, YSAVE, NY2DIM,    &
     		  JAC, WKJAC, NWKJAC, JACPVT, NJCPVT, MONITR,  &
     		  LDERIV, ITASK, ITRACE, IFAIL)
      IF (IFAIL.EQ.1 .OR. IFAIL.GE.14) STOP
      IFAIL = 0
      CALL D02NYF(...)
          .
          .
          .
      STOP
      END

4  References

5  Parameters

1:     NEQ – INTEGERInput

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

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

2:     LDYSAV – INTEGERInput

On entry: a bound on the maximum number of equations to be solved during the integration.

Constraint: $\mathbf{LDYSAV}\ge \mathbf{NEQ}$.

3:     T – REAL (KIND=nag_wp)Input/Output

On 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.

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

4:     TOUT – REAL (KIND=nag_wp)Input/Output

On 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}$.

On exit: normally unchanged. However, when $\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.

5:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput/Output

On 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.

On exit: the computed solution vector, evaluated at T (usually $\mathbf{T}=\mathbf{TOUT}$).

6:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput/Output

On entry: if $\mathbf{LDERIV}\left(1\right)=\mathrm{.TRUE.}$, YDOT must contain approximations to the time derivatives $y^{\prime }$ of the vector $y$.

If $\mathbf{LDERIV}\left(1\right)=\mathrm{.FALSE.}$, YDOT need not be set on entry.

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

7:     RWORK($50+4\times \mathbf{NEQ}$) – REAL (KIND=nag_wp) arrayCommunication Array

8:     RTOL($*$) – REAL (KIND=nag_wp) arrayInput

Note: the dimension of the array RTOL must be at least $1$ if $\mathbf{ITOL}=1$ or $2$, and at least $\mathbf{NEQ}$ otherwise.

On entry: the relative local error tolerance.

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

9:     ATOL($*$) – REAL (KIND=nag_wp) arrayInput

Note: the dimension of the array ATOL must be at least $1$ if $\mathbf{ITOL}=1$ or $3$, and at least $\mathbf{NEQ}$ otherwise.

On entry: the absolute local error tolerance.

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

10:   ITOL – INTEGERInput

On entry: a value to indicate the form of the local error test. ITOL indicates to D02NHF 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 routine.

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

11:   INFORM($23$) – INTEGER arrayCommunication Array

12:   RESID – SUBROUTINE, supplied by the user.External Procedure

RESID must evaluate the residual

$$r=g\left(t,y\right)-A\left(t,y\right)y^{\prime }$$

in one case and

$$r=-A\left(t,y\right)y^{\prime }$$

in another.

The specification of RESID is:

SUBROUTINE RESID ( NEQ, T, Y, YDOT, R, IRES) INTEGER  NEQ, IRES REAL (KIND=nag_wp)  T, Y(NEQ), YDOT(NEQ), R(NEQ)

1:     NEQ – INTEGERInput

On entry: the number of equations being solved.

2:     T – REAL (KIND=nag_wp)Input

On entry: $t$, the current value of the independent variable.

3:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: the value of $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{NEQ}$.

4:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: the value of $y_{\mathit{i}}^{\prime }$, for $\mathit{i}=1,2,\dots ,\mathbf{NEQ}$, at $t$.

5:     R(NEQ) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{R}\left(\mathit{i}\right)$ must contain the $\mathit{i}$th component of $r$, for $\mathit{i}=1,2,\dots ,\mathbf{NEQ}$, where

$$r=g\left(t,y\right)-A\left(t,y\right)y^{\prime }$$

(1)

or

$$r=-A\left(t,y\right)y^{\prime }$$

(2)

and where the definition of $r$ is determined by the input value of IRES.

6:     IRES – INTEGERInput/Output

On entry: the form of the residual that must be returned in array R.

$\mathbf{IRES}=-1$

The residual defined in equation (2) must be returned.

$\mathbf{IRES}=1$

The residual defined in equation (1) must be returned.

On exit: should be unchanged unless one of the following actions is required of the integrator, in which case IRES should be set accordingly.

$\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}$.

$\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$. 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}$.

$\mathbf{IRES}=4$

Indicates to the integrator to stop its current operation and to enter MONITR immediately with parameter $\mathbf{IMON}=-2$.

RESID must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NHF is called. Parameters denoted as Input must not be changed by this procedure.

13:   YSAV(LDYSAV,SDYSAV) – REAL (KIND=nag_wp) arrayCommunication Array

14:   SDYSAV – INTEGERInput

On entry: an appropriate value for SDYSAV is described in the specifications of the integrator setup routines D02MVF, D02NVF and D02NWF. This value must be the same as that supplied to the integrator setup routine.

15:   JAC – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

JAC must evaluate the Jacobian of the system. If this option is not required, the actual argument for JAC must be the dummy routine D02NHZ. (D02NHZ is included in the NAG Library.) You must indicate to the integrator whether this option is to be used by setting the parameter JCEVAL appropriately in a call to the banded linear algebra setup routine D02NTF.

First we must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, $y^{\prime }$, generated internally, has the form

$$y^{\prime }=\left(y-z\right)/\left(hd\right)\text{,}$$

where $h$ is the current step size and $d$ is a parameter that depends on the integration method in use. The vector $y$ is the current solution and the vector $z$ depends on information from previous time steps. This means that $\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\left(t,y\right)y^{\prime }-g\left(t,y\right)=0$$

but it is solved in the form

$$r\left(t,y\right)=0\text{,}$$

where $r$ is the function defined by

$$r\left(t,y\right)=\left(hd\right)\left(A\left(t,y\right)\left(y-z\right)/\left(hd\right)-g\left(t,y\right)\right)\text{.}$$

It is the Jacobian matrix $\frac{\partial r}{\partial y}$ that you must supply in JAC as follows:

$$\frac{\partial r_i}{\partial y_j}=a_{ij}\left(t,y\right)+\left(hd\right)\frac{\partial }{\partial y_j}\left(\sum _{k=1}^{\mathbf{NEQ}}{a}_{ik}\left(t,y\right)y_k^{\prime }-g_i\left(t,y\right)\right)\text{.}$$

The specification of JAC is:

SUBROUTINE JAC ( NEQ, T, Y, YDOT, H, D, ML, MU, P) INTEGER  NEQ, ML, MU REAL (KIND=nag_wp)  T, Y(NEQ), YDOT(NEQ), H, D, P(ML+MU+1,NEQ)

1:     NEQ – INTEGERInput

On entry: the number of equations being solved.

2:     T – REAL (KIND=nag_wp)Input

On entry: $t$, the current value of the independent variable.

3:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{NEQ}$, the current solution component.

4:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: the derivative of the solution at the current point $t$.

5:     H – REAL (KIND=nag_wp)Input

On entry: the current step size.

6:     D – REAL (KIND=nag_wp)Input

On entry: the parameter $d$ which depends on the integration method.

7:     ML – INTEGERInput

8:     MU – INTEGERInput

On entry: the number of subdiagonals and superdiagonals respectively in the band.

9:     P($\mathbf{ML}+\mathbf{MU}+1$,NEQ) – REAL (KIND=nag_wp) arrayInput/Output

On entry: is set to zero.

On exit: elements of the Jacobian matrix $\frac{\partial r}{\partial y}$ stored as specified by the following pseudocode

 DO 20 I = 1, NEQ J1 = MAX(I-ML,1) J2 = MIN(I+MU,NEQ) DO 10 J = J1, J2 K = MIN(ML+1-I,0)+J P(K,I) = δR/δY(I,J) 10 CONTINUE 20 CONTINUE 

See also the routine document for F07BDF (DGBTRF).

Only nonzero elements of this array need be set, since it is preset to zero before the call to JAC.

JAC must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NHF is called. Parameters denoted as Input must not be changed by this procedure.

16:   WKJAC(NWKJAC) – REAL (KIND=nag_wp) arrayCommunication Array

17:   NWKJAC – INTEGERInput

On entry: the dimension of the array WKJAC as declared in the (sub)program from which D02NHF is called. This value must be the same as that supplied to the linear algebra setup routine D02NTF.

Constraint: $\mathbf{NWKJAC}\ge \left(2m_L+m_U+1\right)\times \mathbf{NEQ}$, where $m_L$ and $m_U$ are the number of subdiagonals and superdiagonals respectively in the band, defined by a call to D02NTF.

18:   JACPVT(NJCPVT) – INTEGER arrayCommunication Array

19:   NJCPVT – INTEGERInput

On entry: the size of the array JACPVT. This value must be the same as that supplied to the linear algebra setup routine D02NTF.

Constraint: $\mathbf{NJCPVT}\ge \mathbf{NEQ}$.

20:   MONITR – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

MONITR performs tasks requested by you. If this option is not required, then the actual argument for MONITR must be the dummy routine D02NBY. (D02NBY is included in the NAG Library.)

The specification of MONITR is:

SUBROUTINE MONITR ( NEQ, LDYSAV, T, HLAST, HNEXT, Y, YDOT, YSAV, R, ACOR, IMON, INLN, HMIN, HMAX, NQU) INTEGER  NEQ, LDYSAV, IMON, INLN, NQU REAL (KIND=nag_wp)  T, HLAST, HNEXT, Y(NEQ), YDOT(NEQ), YSAV(LDYSAV,$\mathit{sdysav}$), R(NEQ), ACOR(NEQ,2), HMIN, HMAX

where $\mathit{sdysav}$ is the numerical value of SDYSAV in the call of D02NHF.

1:     NEQ – INTEGERInput

On entry: the number of equations being solved.

2:     LDYSAV – INTEGERInput

On entry: an upper bound on the number of equations to be solved.

3:     T – REAL (KIND=nag_wp)Input

On entry: the current value of the independent variable.

4:     HLAST – REAL (KIND=nag_wp)Input

On entry: the last step size successfully used by the integrator.

5:     HNEXT – REAL (KIND=nag_wp)Input/Output

On entry: the step size that the integrator proposes to take on the next step.

On exit: the next step size to be used. If this is different from the input value, then IMON must be set to $4$.

6:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput/Output

On entry: $y$, the values of the dependent variables evaluated at $t$.

On exit: these values must not be changed unless IMON is set to $2$.

7:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: the time derivatives $y^{\prime }$ of the vector $y$.

8:     YSAV(LDYSAV,$\mathit{sdysav}$) – REAL (KIND=nag_wp) arrayInput

On entry: workspace to enable you to carry out interpolation using either of the routines D02XJF or D02XKF.

9:     R(NEQ) – REAL (KIND=nag_wp) arrayInput

On entry: if $\mathbf{IMON}=0$ and $\mathbf{INLN}=3$, then the first NEQ elements contain the residual vector $A\left(t,y\right)y^{\prime }-g\left(t,y\right)$.

10:   ACOR(NEQ,$2$) – REAL (KIND=nag_wp) arrayInput

On entry: with $\mathbf{IMON}=1$, $\mathbf{ACOR}\left(i,1\right)$ contains the weight used for the $i$th equation when the norm is evaluated, and $\mathbf{ACOR}\left(i,2\right)$ contains the estimated local error for the $i$th equation. The scaled local error at the end of a timestep may be obtained by calling the real function D02ZAF as follows:

 IFAIL = 1 ERRLOC = D02ZAF(NEQ, ACOR(1,2), ACOR(1,1), IFAIL) ! CHECK IFAIL BEFORE PROCEEDING 

11:   IMON – INTEGERInput/Output

On entry: a flag indicating under what circumstances MONITR was called:

$\mathbf{IMON}=-2$

Entry from the integrator after $\mathbf{IRES}=4$ (set in RESID) caused an early termination (this facility could be used to locate discontinuities).

$\mathbf{IMON}=-1$

The current step failed repeatedly.

$\mathbf{IMON}=0$

Entry after a call to the internal nonlinear equation solver (see INLN).

$\mathbf{IMON}=1$

The current step was successful.

On exit: may be reset to determine subsequent action in D02NHF.

$\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}$.

$\mathbf{IMON}=-1$

Allow the integrator to continue with its own internal strategy. The integrator will try up to three restarts unless $\mathbf{IMON}\ne -1$ on exit.

$\mathbf{IMON}=0$

Return to the internal nonlinear equation solver, where the action taken is determined by the value of INLN (see INLN).

$\mathbf{IMON}=1$

Normal exit to the integrator 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 MONITR, 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 the call to MONITR. 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.

12:   INLN – INTEGEROutput

On exit: the action to be taken by the internal nonlinear equation solver when MONITR is exited with $\mathbf{IMON}=0$. By setting $\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.

13:   HMIN – REAL (KIND=nag_wp)Input/Output

On entry: the minimum step size to be taken on the next step.

On exit: the minimum step size to be used. If this is different from the input value, then IMON must be set to $3$ or $4$.

14:   HMAX – REAL (KIND=nag_wp)Input/Output

On entry: the maximum step size to be taken on the next step.

On exit: the maximum step size to be used. If this is different from the input value, then IMON must be set to $3$ or $4$. If HMAX is set to zero, no limit is assumed.

15:   NQU – INTEGERInput

On entry: the order of the integrator used on the last step. This is supplied to enable you to carry out interpolation using either of the routines D02XJF or D02XKF.

MONITR must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NHF is called. Parameters denoted as Input must not be changed by this procedure.

21:   LDERIV($2$) – LOGICAL arrayInput/Output

On entry: $\mathbf{LDERIV}\left(1\right)$ must be set to .TRUE. if you have supplied both an initial $y$ and an initial $y^{\prime }$. $\mathbf{LDERIV}\left(1\right)$ must be set to .FALSE. if only the initial $y$ has been supplied.

$\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$ and $y^{\prime }$. Note that $y$ and $y^{\prime }$, if supplied, are used as initial estimates. This method involves taking a small step at the start of the integration, and if $\mathbf{ITASK}=6$ on entry, T and TOUT will be set to the result of taking this small step. $\mathbf{LDERIV}\left(2\right)$ must be set to .FALSE. if the integrator is to use functional iteration to evaluate the initial $y$ and $y^{\prime }$, and if this fails a modified Newton method will then be attempted. $\mathbf{LDERIV}\left(2\right)=\mathrm{.TRUE.}$ is recommended if there are implicit equations or the initial $y$ and $y^{\prime }$ are zero.

On exit: $\mathbf{LDERIV}\left(1\right)$ is normally unchanged. However if $\mathbf{ITASK}=6$ and internal initialization was successful then $\mathbf{LDERIV}\left(1\right)=\mathrm{.TRUE.}$.

$\mathbf{LDERIV}\left(2\right)=\mathrm{.TRUE.}$, if implicit equations were detected. Otherwise $\mathbf{LDERIV}\left(2\right)=\mathrm{.FALSE.}$.

22:   ITASK – INTEGERInput

On 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 routines before the first call to the integrator, or specified in the optional input routine before a continuation call. TCRIT 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. TCRIT must be specified as under $\mathbf{ITASK}=4$.

$\mathbf{ITASK}=6$

The integrator will solve for the initial values of $y$ and $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$ and $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$ will have been advanced a short distance from the initial point.

Note:  if D02NHF is recalled with a different value of ITASK (and TOUT altered), then the initialization procedure is repeated, possibly leading to different initial conditions.

Constraint: $1\le \mathbf{ITASK}\le 6$.

23:   ITRACE – INTEGERInput

On 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 X04AAF).

$\mathbf{ITRACE}>0$

Warning messages are printed as above, and on the current advisory message unit (see X04ABF) 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.

24:   IFAIL – INTEGERInput/Output

On entry: IFAIL must be set to $0$, $-1\text{ or }1$. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{ or }1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if $\mathbf{IFAIL}\ne \mathbf{0}$ on exit, the recommended value is $-1$. When the value $-\mathbf{1}\text{ or }1$ is used it is essential to test the value of IFAIL on exit.

On exit: $\mathbf{IFAIL}=\mathbf{0}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

$\mathbf{IFAIL}=1$

An illegal input was detected on entry, or after an internal call to MONITR. If $\mathbf{ITRACE}>-1$, then the form of the error will be detailed on the current error message unit (see X04AAF).

$\mathbf{IFAIL}=2$

The maximum number of steps specified has been taken (see the description of optional inputs in the integrator setup routines and the optional input continuation routine, D02NZF).

$\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.

$\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.

$\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.

$\mathbf{IFAIL}=6$

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

$\mathbf{IFAIL}=7$

RESID set its error flag ($\mathbf{IRES}=3$) continually despite repeated attempts by the integrator to avoid this.

$\mathbf{IFAIL}=8$

$\mathbf{LDERIV}\left(1\right)=\mathrm{.FALSE.}$ on entry but the internal initialization routine was unable to initialize $y^{\prime }$ (more detailed information may be directed to the current error message unit, see X04AAF).

$\mathbf{IFAIL}=9$

A singular Jacobian $\frac{\partial r}{\partial y}$ has been encountered. 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 X04AAF).

$\mathbf{IFAIL}=11$

RESID signalled the integrator to halt the integration and return ($\mathbf{IRES}=2$). Integration was successful as far as T.

$\mathbf{IFAIL}=12$

MONITR set $\mathbf{IMON}=-2$ and so forced a return but the integration was successful as far as T.

$\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 D02NHF is unable to start the integration.

$\mathbf{IFAIL}=15$

The linear algebra setup routine D02NTF was not called before the call to D02NHF.

7  Accuracy

8  Further Comments

9  Example

9.1  Program Text

Program Text (d02nhfe.f90)

9.2  Program Data

Program Data (d02nhfe.d)

9.3  Program Results

Program Results (d02nhfe.r)