NAG Library Chapter Introduction

D02 (ode)
Ordinary Differential Equations

1
Scope of the Chapter

This chapter is concerned with the numerical solution of ordinary differential equations. There are two main types of problem: those in which all boundary conditions are specified at one point (initial value problems), and those in which the boundary conditions are distributed between two or more points (boundary value problems and eigenvalue problems). Routines are available for initial value problems, two-point boundary value problems and Sturm–Liouville eigenvalue problems.

2
Background to the Problems

For most of the routines in this chapter a system of ordinary differential equations must be written in the form
y1=f1x,y1,y2,,yn, y2=f2x,y1,y2,,yn, yn=fnx,y1,y2,,yn,  
that is the system must be given in first-order form. The n dependent variables (also, the solution) y1,y2,,yn are functions of the independent variable x, and the differential equations give expressions for the first derivatives yi=dyi dx  in terms of x and y1,y2,,yn. For a system of n first-order equations, n associated boundary conditions are usually required to define the solution.
A more general system may contain derivatives of higher order, but such systems can almost always be reduced to the first-order form by introducing new variables. For example, suppose we have the third-order equation
z+zz+kl-z2=0.  
We write y1=z, y2=z, y3=z, and the third-order equation may then be written as the system of first-order equations
y1=y2 y2=y3 y3=-y1y3-kl-y22.  
For this system n=3 and we require 3 boundary conditions in order to define the solution. These conditions must specify values of the dependent variables at certain points. For example, we have an initial value problem if the conditions are
y1=0.1  at ​x=0 y2=0.1  at ​x=0 y3=0.1  at ​x=0.  
These conditions would enable us to integrate the equations numerically from the point x=0 to some specified end point. We have a boundary value problem if the conditions are
y1=0  at ​x=0 y2=0  at ​x=0 y2=1  at ​x=10.  
These conditions would be sufficient to define a solution in the range 0x10, but the problem could not be solved by direct integration (see Section 2.2). More general boundary conditions are permitted in the boundary value case.
It is sometimes advantageous to solve higher-order systems directly. In particular, there is an initial value routine to solve a system of second-order ordinary differential equations of the special form
y1=f1x,y1,y2,,yn, y2=f2x,y1,y2,,yn, yn=fnx,y1,y2,,yn.  
For this second-order system initial values of the derivatives of the dependent variables, yi, for i=1,2,,n, are required.
There is also a boundary value routine that can treat directly a mixed order system of ordinary differential equations.
There is a broader class of initial value problems known as differential algebraic systems which can be treated. Such a system may be defined as
y = fx,y,z 0 = gx,y,z  
where y and f are vectors of length n and g and z are vectors of length m. The functions g represent the algebraic part of the system.
In addition implicit systems can also be solved, that is systems of the form
Ax,yy=fx,y  
where A is a matrix of functions; such a definition can also incorporate algebraic equations. Note that general systems of this form may contain higher-order derivatives and that they can usually be transformed to first-order form, as above.

2.1
Initial Value Problems

To solve first-order systems, initial values of the dependent variables yi, for i=1,2,,n, must be supplied at a given point, a. Also a point, b, at which the values of the dependent variables are required, must be specified. The numerical solution is then obtained by a step-by-step calculation which approximates values of the variables yi, for i=1,2,,n, at finite intervals over the required range a,b. The routines in this chapter adjust the step length automatically to meet specified accuracy tolerances. Although the accuracy tests used are reliable over each step individually, in general an accuracy requirement cannot be guaranteed over a long range. For many problems there may be no serious accumulation of error, but for unstable systems small perturbations of the solution will often lead to rapid divergence of the calculated values from the true values. A simple check for stability is to carry out trial calculations with different tolerances; if the results differ appreciably the system is probably unstable. Over a short range, the difficulty may possibly be overcome by taking sufficiently small tolerances, but over a long range it may be better to try to reformulate the problem.
A special class of initial value problems are those for which the solutions contain rapidly decaying transient terms. Such problems are called stiff; an alternative way of describing them is to say that certain eigenvalues of the Jacobian matrix fi yj  have large negative real parts when compared to others. These problems require special methods for efficient numerical solution; the methods designed for non-stiff problems when applied to stiff problems tend to be very slow, because they need small step lengths to avoid numerical instability. A full discussion is given in Hall and Watt (1976) and a discussion of the methods for stiff problems is given in Berzins et al. (1988).

2.2
Boundary Value Problems

In general, a system of nonlinear differential equations with boundary conditions at two or more points cannot be guaranteed to have a solution. The solution, if it exists, has to be determined iteratively. A comprehensive treatment of the numerical solution of boundary value problems can be found in Ascher et al. (1988) and Keller (1992). The methods for this chapter are discussed in Ascher et al. (1979), Ascher and Bader (1987) and Gladwell (1987).

2.2.1
Collocation methods

In the collocation method, the solution components are approximated by piecewise polynomials on a mesh. The coefficients of the polynomials form the unknowns to be computed. The approximation to the solution must satisfy the boundary conditions and the differential equations at collocation points in each mesh sub-interval. A modified Newton method is used to solve the nonlinear equations. The mesh is refined by trying to equidistribute the estimated error over the whole interval. An initial estimate of the solution across the mesh is required.

2.2.2
Shooting methods

In the shooting method, the unknown boundary values at the initial point are estimated to form an initial value problem, and the equations are then integrated to the final point. At the final point the computed solution and the known boundary conditions should be equal. The condition for equality gives a set of nonlinear equations for the estimated values, which can be solved by Newton's method or one of its variants. The iteration cannot be guaranteed to converge, but it is usually successful if
It may be necessary to simplify the problem and carry out some preliminary calculations, in order to obtain suitable starting values. A fuller discussion is given in Chapters 16, 17 and 18 of Hall and Watt (1976), Chapter 11 of Gladwell and Sayers (1980) and Chapter 8 of Gladwell (1979a).

2.2.3
Finite difference methods

If a boundary value problem seems insoluble by the above method and a good estimate for the solution of the problem is known at all points of the range then a finite difference method may be used. Finite difference equations are set up on a mesh of points and estimated values for the solution at the grid points are chosen. Using these estimated values as starting values a Newton iteration is used to solve the finite difference equations. The accuracy of the solution is then improved by deferred corrections or the addition of points to the mesh or a combination of both. The method does not suffer from the difficulties associated with the shooting method but good initial estimates of the solution may be required in some cases and the method is unlikely to be successful when the solution varies very rapidly over short ranges. A discussion is given in Chapters 9 and 11 of Gladwell and Sayers (1980) and Chapter 4 of Gladwell (1979a).

2.3
Chebyshev Collocation for Linear Differential Equations

The collocation method gives a different approach to the solution of ordinary differential equations. It can be applied to problems of either initial value or boundary value type. Suppose the approximate solution is represented in polynomial form, say as a series of Chebyshev polynomials. The coefficients may be determined by matching the series to the boundary conditions, and making it satisfy the differential equation at a number of selected points in the range. The calculation is straightforward for linear differential equations (nonlinear equations may also be solved by an iterative technique based on linearization). The result is a set of Chebyshev coefficients, from which the solution may be evaluated at any point using e02akf. A fuller discussion is given in Chapter 24 of Gladwell (1979a) and Chapter 11 of Gladwell and Sayers (1980).
This method can be useful for obtaining approximations to standard mathematical functions. For example, suppose we require values of the Bessel function J13x over the range 0,5, for use in another calculation. We solve the Bessel differential equation by collocation and obtain the Chebyshev coefficients of the solution, which we can use to construct a function for J13x. (Note that routines for many common standard functions are already available in Chapter S.)

2.4
Eigenvalue Problems

Sturm–Liouville problems of the form
px y + q x,λ y=0  
with appropriate boundary conditions given at two points, can be solved by a Scaled Prüfer method. In this method the differential equation is transformed to another which can be solved for a specified eigenvalue by a shooting method. A discussion is given in Chapter 11 of Gladwell and Sayers (1980) and a complete description is given in Pryce (1986). Some more general eigenvalue problems can be solved by the methods described in Section 2.2.

3
Recommendations on Choice and Use of Available Routines

There are no routines which deal directly with complex equations. These may however be transformed to larger systems of real equations of the required form. Split each equation into its real and imaginary parts and solve for the real and imaginary parts of each component of the solution. Whilst this process doubles the size of the system and may not always be appropriate it does make available for use the full range of routines provided presently.

3.1
Initial Value Problems

In general, for non-stiff first-order systems, Runge–Kutta (RK) routines should be used. For the usual requirement of integrating across a range the appropriate routines are d02pef and d02pqf; d02pqf is a setup routine for d02pef. For more complex tasks there are forward and reverse communication variants (Section 3.3.3 in How to Use the NAG Library and its Documentation) of single step routines with corresponding interpolator; for direct communication these are d02pff and d02psf, while for reverse communication these are d02pgf, d02phf and d02pjf. There are also related utility routines d02prf, d02ptf and d02puf. When a system is to be integrated over a long range or with relatively high accuracy requirements the variable-order, variable-step Adams' codes may be more efficient. The appropriate routine in this case is d02cjf. For more complex tasks using an Adams' code there are a further six related routines: d02qff, d02qgf, d02qwf, d02qxf, d02qyf and d02qzf.
For stiff systems, that is those which usually contain rapidly decaying transient components, the Backward Differentiation Formula (BDF) variable-order, variable-step codes should be used. The appropriate routine in this case is d02ejf. For more complex tasks where the system residual is difficult to evaluate in direct communication, or is coupled with algebraic equations, there are a collection of routines in Sub-chapter D02M–N. These routines can treat implicit differential algebraic systems, they also contain additional methods (beyond BDF techniques) which may be appropriate in some circumstances.
If you are not sure how to classify a problem, you are advised to perform some preliminary calculations with d02pef, which can indicate whether the system is stiff. We also advise performing some trial calculations with d02pef (RK), d02cjf (Adams) and d02ejf (BDF) so as to determine which type of routine is best applied to the problem. The conclusions should be based on the computer time used and the number of evaluations of the derivative function fi. See Gladwell (1979b) for more details.
For second-order systems of the special form described in Section 2 the Runge–Kutta–Nystrom (RKN) routine d02laf should be used.

3.1.1
Runge–Kutta routines

The basic RK routines are d02pff (direct communication) and d02pgf (reverse communication) which take one integration step at a time. An alternative to d02pff is d02pef, which provides output at user-specified points. The initialization of d02pef, d02pff or d02pgf and the setting of optional inputs, including choice of method, is made by a call to the setup routine d02pqf. Optional output information about the integration and about error assessment, if selected, can be obtained by calls to the diagnostic routines d02ptf and d02puf respectively. d02psf may be used to interpolate on information produced by d02pff to give solution and derivative values between the integration points. Similarly d02phf may be used to setup an interpolator on information produced by d02pgf, and d02pjf can evaluate that interpolator to give solution and derivative values between integration points; these routines are recommended when a high-order RK method is specified in the setup routine. d02prf may be used to reset the end of the integration range whilst integrating using d02pff or d02pgf.
There is a simple driving routine d02bjf, which integrates a system over a range and, optionally, computes intermediate output and/or determines the position where a specified function of the solution is zero.
For well-behaved systems with no difficulties such as stiffness or singularities, the Merson form of the RK method, as used by d02bgf, works well for low accuracy calculations. d02bhf also uses the Merson form and can additionally find a root of a supplied equation involving solution components.

3.1.2
Adams' routines

The general Adams' variable-order variable-step routine is d02qff, which provides a choice of automatic error control and the option of a sophisticated root-finding technique. Reverse communication for both the differential equation and root definition function is provided in d02qgf, which otherwise has the same facilities as d02qff. A reverse communication routine makes a return to the calling subroutine for evaluations of equations rather than calling a user-supplied procedure. The initialization of either of d02qff and d02qgf and the setting of optional inputs is made by a call to the setup routine d02qwf. Optional output information about the integration and any roots detected can be obtained by calls to the diagnostic routines d02qxf and d02qyf respectively. d02qzf may be used to interpolate on information produced by d02qff or d02qgf to give solution and derivative values between the integration points.
There is a simple driving routine d02cjf, which integrates a system over a range and, optionally, computes intermediate output and/or determines the position where a specified function of the solution is zero.

3.1.3
BDF routines

General routines for explicit and implicit ordinary differential equations with a wide range of options for integrator choice and special forms of numerical linear algebra are provided in Sub-chapter D02M–N. A separate document describing the use of this sub-chapter is given immediately before the routines of the sub-chapter.
There are three utility routines available for use with Sub-chapter D02M–N routines. d02xjf and d02xkf can be used to interpolate to a solution at a given point using the natural and C1 interpolants respectively. d02zaf can be used to return the weighted norm of the local error estimate calculated by Sub-chapter D02M–N routines.
There is a simple driving routine d02ejf, which integrates a system over a range and, optionally, computes intermediate output and/or determines the position where a specified function of the solution is zero. It has a specification similar to the Adams' routine d02cjf except that to solve the equations arising in the BDF method an approximation to the Jacobian fi yj  is required. This approximation can be calculated internally but you may supply an analytic expression. In most cases supplying a correct analytic expression will reduce the amount of computer time used.

3.1.4
Runge–Kutta–Nystrom routines

The Runge–Kutta–Nystrom routine d02laf uses either a low- or high-order method (chosen by you). The choice of method and error control and the setting of optional inputs is made by a call to the setup routine d02lxf. Optional output information about the integration can be obtained by a call to the diagnostic routine d02lyf. When the low-order method has been employed d02lzf may be used to interpolate on information produced by d02laf to give the solution and derivative values between the integration points.

3.2
Boundary Value Problems

In general, for a nonlinear system of mixed order with separated boundary conditions, the collocation method (d02tlf and its associated routines) can be used. Problems of a more general nature can often be transformed into a suitable form for treatment by d02tlf, for example nonseparated boundary conditions or problems with unknown parameters (see Section 9 in d02tvf for details).
For simple boundary value problems with assigned boundary values you may prefer to use a code based on the shooting method or finite difference method for which there are routines with simple calling sequences (d02haf and d02gaf).
For difficult boundary value problems, where you need to exercise some control over the calculation, and where the collocation method proves unsuccessful, you may wish to try the alternative methods of shooting (d02saf) or finite differences (d02raf).
Note that it is not possible to make a fully automatic boundary value routine, and you should be prepared to experiment with different starting values or a different routine if the problem is at all difficult.

3.2.1
Collocation methods

The collocation routine d02tlf solves a nonlinear system of mixed order boundary value problems with separated boundary conditions. The initial mesh and accuracy requirements must be specified by a call to the setup routine d02tvf. Optional output information about the final mesh and estimated maximum error can be obtained by a call to the diagnostic routine d02tzf. The solution anywhere on the mesh can be computed by a call to the interpolation routine d02tyf. If d02tlf is being used to solve a sequence of related problems then the continuation routine d02txf should also be used.

3.2.2
Shooting methods

d02haf may be used for simple boundary value problems, where the unknown arguments are the missing boundary conditions. More general boundary value problems may be handled by using d02hbf. This routine allows for a generalized argument structure, and is fairly complicated. The older routine d02agf has been retained for use when an interior matching-point is essential; otherwise the newer routine d02hbf should be preferred.
For particularly complicated problems where, for example, the arguments must be constrained or the range of integration must be split to enable the shooting method to succeed, the recommended routine is d02saf, which extends the facilities provided by d02hbf. If you are an experienced user d02saf permits you much more control over the calculation than does d02hbf; in particular you are permitted precise control of solution output and intermediate monitoring information.

3.2.3
Finite difference methods

d02gaf may be used for simple boundary value problems with assigned boundary values. The calling sequence of d02gaf is very similar to that for d02haf discussed above.
You may find that convergence is difficult to achieve using d02gaf since only specifying the unknown boundary values and the position of the finite difference mesh is permitted. In such cases you may use d02raf, which permits specification of an initial estimate for the solution at all mesh points and allows the calculation to be influenced in other ways too. d02raf is designed to solve a general nonlinear two-point boundary value problem with nonlinear boundary conditions.
A routine, d02gbf, is also supplied specifically for the general linear two-point boundary value problem written in a standard ‘textbook’ form.
You are advised to use interpolation routines from Chapter E01 to obtain solution values at points not on the final mesh.

3.2.4
Chebyshev integration method

The Chebyshev integration method is an implementation of the Chebyshev collocation method (see Section 3.3) which is fully described and compared against other implementations in Muite (2010). d02uef solves a linear constant coefficient boundary value problem using the Chebyshev integration formulation on a Chebyshev Gauss–Lobatto grid and solving in the coefficient space. The required Chebyshev Gauss–Lobatto grid points on a given arbitrary interval a,b can first be generated using d02ucf. Then d02uaf obtains the Chebyshev coefficients for the right-hand side (of system) function discretized on the obtained Chebyshev Gauss–Lobatto grid. d02uef then solves the problem in Chebyshev coefficient space using the integration formulation. Finally d02ubf evaluates the solution (or one of its lower order derivatives) from the set of Chebyshev coefficients returned by d02uef on the Chebyshev Gauss–Lobatto grid on a,b. The set of routines can be used to solve up to fourth order boundary value problems.

3.3
Chebyshev Collocation Method

d02tgf may be used to obtain the approximate solution of a system of differential equations in the form of a Chebyshev series. The routine treats linear differential equations directly, and makes no distinction between initial value and boundary value problems. This routine is appropriate for problems where it is known that the solution is smooth and well-behaved over the range, so that each component can be represented by a single polynomial. Singular problems can be solved using d02tgf as long as their polynomial-like solutions are required.
d02tgf permits the differential equations to be specified in higher order form; that is without conversion to a first-order system. This type of specification leads to a complicated calling sequence. If you are an inexperienced user two simpler routines are supplied. d02jaf solves a single regular linear differential equation of any order whereas d02jbf solves a system of regular linear first-order differential equations.

3.4
Eigenvalue Problems

Two routines, d02kaf and d02kdf, may be used to find the eigenvalues of second-order Sturm–Liouville problems. d02kaf is designed to solve simple problems with regular boundary conditions. d02kaf calls d02kdf, which is designed to solve more difficult problems, for example with singular boundary conditions or on infinite ranges or with discontinuous coefficients.
If the eigenfunctions of the Sturm–Liouville problem are also required, d02kef should be used. (d02kef solves the same types of problem as d02kdf.)

3.5
Summary of Recommended Routines

Problem Routine
  RK Method Adams' Method BDF Method
Initial Value Problems
Driver Routines
     
Integration over a range with optional intermediate output and optional determination of position where a function of the solution becomes zero d02bjf d02cjf d02ejf
Integration of a range with intermediate output d02bjf d02cjf d02ejf
Integration of a range until function of solution becomes zero d02bjf d02cjf d02ejf
Comprehensive Integration Routines d02pef, d02pff, d02pgf, d02phf, d02pjf, d02pqf, d02prf, d02psf, d02ptf and d02puf d02qff, d02qgf, d02qwf, d02qxf, d02qyf and d02qzf Sub-chapter D02M–N routines, d02xjf, d02xkf and d02zaf
Package for Solving Stiff Equations Sub-chapter D02M–N routines
Package for Solving Second-order Systems of Special Form D02L routines
Boundary Value Problems
Collocation Method, Mixed Order
d02tlf, d02tvf, d02txf, d02tyf and d02tzf
Boundary Value Problems
Shooting Method
 
simple argument d02haf
generalized arguments d02agf and d02hbf
additional facilities d02saf
Boundary Value Problems
Finite Difference Method
 
simple argument d02gaf
linear problem d02gbf
full nonlinear problem d02raf
Boundary Value Problems
Chebyshev Collocation, Integration Formulation
 
single linear equation d02uef with d02uaf, d02ubf, d02ucf
Chebyshev Collocation, Linear Problems  
single equation d02jaf
first-order system d02jbf
general system d02tgf
Sturm–Liouville Eigenvalue Problems  
regular problems d02kaf
general problems d02kdf
eigenfunction calculation d02kef

4
Decision Trees

Tree 1: Initial Value Problems

Is the problem first order?   Is the problem known to be stiff?   Use routines described in the D02M–N Sub-chapter Introduction, or their simple driver d02ejf.
yesyes
  no   no
Backward Differentiation Formula: Use routines described in the D02M–N Sub-chapter Introduction, or their simple driver d02ejf.
Adams' method with driver routine: d02cjf
Adams' method with comprehensive suite: d02qff, d02qgf, d02qwf, d02qxf, d02qyf and d02qzf
Runge–Kutta method with driver routine: d02bjf
Runge–Kutta method with comprehensive suite: d02pef, d02pff, d02pgf, d02phf, d02pjf, d02pqf, d02prf, d02psf, d02ptf and d02puf
Is the problem of the form: y=fx,y   Use the D02L routines
yes
  no
Convert to first order problem: is the problem known to be stiff?   Use routines described in the D02M–N Sub-chapter Introduction, or their simple driver d02ejf.
yes
  no
Backward Differentiation Formula: use routines described in the D02M–N Sub-chapter Introduction, or their simple driver d02ejf.
Adams' method with driver routine: d02cjf
Adams' method with comprehensive suite: d02qff, d02qgf, d02qwf, d02qxf, d02qyf and d02qzf
Runge–Kutta method with driver routine: d02bjf
Runge–Kutta method with comprehensive suite: d02pef, d02pff, d02pgf, d02phf, d02pjf, d02pqf, d02prf, d02psf, d02ptf and d02puf

Tree 2: Boundary Value Problems

Is the problem simply of the form y'=fx,y?   Are only boundary values to be determined?   Shooting method: d02haf
Finite differences: d02gaf
Collocation: d02tlf
yesyes
  no   no
Shooting method: d02hbf
Finite differences: d02gbf
Collocation, piecewise polynomials: d02tlf
Collocation, Chebyshev polynomials: d02jaf, d02jbf or d02uef
Shooting method: d02saf
Finite differences: d02raf
Collocation, piecewise polynomials: d02tlf
Collocation, Chebyshev polynomials: d02jaf, d02jbf or d02uef

5
Functionality Index

Differentiation of a function discretized on Chebyshev Gauss–Lobatto points d02udf
Linear constant coefficient boundary value problem, 
    Chebyshev spectral integration method, 
        Chebyshev coefficients generator for a function discretized on Chebyshev Gauss–Lobatto grid d02uaf
        Chebyshev coefficients to function values on Chebyshev Gauss–Lobatto grid d02ubf
        Chebyshev Gauss–Lobatto grid generator d02ucf
        Chebyshev integration solver for linear constant coefficient boundary value problem d02uef
        Clenshaw–Curtis quadrature weights generator at Chebyshev Gauss–Lobatto points d02uyf
        Evaluation on uniform grid of function by Barycentric Lagrange interpolation d02uwf
        value of kth Chebyshev polynomial d02uzf
Second-order Sturm–Liouville problems, 
    regular/singular system, finite/infinite range, 
        eigenvalue and eigenfunction d02kef
        eigenvalue only d02kdf
    regular system, finite range, user-specified break-points, 
        eigenvalue only d02kaf
System of first-order ordinary differential equations, initial value problems, 
    C1-interpolant d02xkf
    comprehensive integrator routines for stiff systems, 
        continuation to call d02nef d02mcf
        explicit ordinary differential equations, 
            banded Jacobian d02ncf
            full Jacobian d02nbf
            sparse Jacobian d02ndf
        explicit ordinary differential equations (reverse communication): 
            full Jacobian d02nmf
        implicit ordinary differential equations coupled with algebraic equations, 
            banded Jacobian d02nhf
            banded Jacobian selector for DASSL integrator d02npf
            DASSL integrator d02nef
            full Jacobian d02ngf
            integrator setup for DASSL d02mwf
            sparse Jacobian d02njf
        implicit ordinary differential equations coupled with algebraic equations (reverse communication) d02nnf
    comprehensive integrator routines using Adams' method with root-finding option, 
        diagnostic routine d02qxf
        diagnostic routine for root-finding d02qyf
        direct communication d02qff
        interpolant d02qzf
        reverse communication d02qgf
        setup routine d02qwf
    comprehensive integrator routines using Runge–Kutta methods, 
        diagnostic routine d02ptf
        diagnostic routine for global error assessment d02puf
        interpolant, reverse communication d02phf
        interpolant and interpolation, direct communication d02psf
        interpolation, reverse communication d02pjf
        over a range with intermediate output d02pef
        over a step, direct communication d02pff
        over a step, reverse communication d02pgf
        reset end of range d02prf
        setup routine d02pqf
    compute weighted norm of local error estimate d02zaf
    enquiry routine for use with sparse Jacobian d02nrf
    integrator diagnostic routine d02nyf
    integrator setup for backward differentiation formulae method for SPRINT integrator d02nvf
    integrator setup for Blend method for SPRINT integrator d02nwf
    integrator setup for DASSL method for SPRINT integrator d02mvf
    linear algebra diagnostic routine for sparse Jacobians d02nxf
    linear algebra setup for banded Jacobians d02ntf
    linear algebra setup for full Jacobians d02nsf
    linear algebra setup for sparse Jacobians d02nuf
    natural interpolant d02mzf
    natural interpolant (for use by MONITR subroutine) d02xjf
    setup routine for continuation calls to integrator d02nzf
    simple driver routines, 
        Runge–Kutta–Merson method, 
            until a function of the solution is zero d02bhf
            until a specified component attains a given value d02bgf
        Runge–Kutta method, 
            until (optionally) a function of the solution is zero, with optional intermediate output d02bjf
        variable-order variable-step Adams' method, 
            until (optionally) a function of the solution is zero, with optional intermediate output d02cjf
        variable-order variable-step backward differentiation formulae method for stiff systems, 
            until (optionally) a function of the solution is zero, with optional intermediate output d02ejf
System of ordinary differential equations, boundary value problems, 
    collocation and least squares, 
        single nth-order linear equation d02jaf
        system of first-order linear equations d02jbf
        system of nth-order linear equations d02tgf
    comprehensive routines using a collocation technique, 
        continuation routine d02txf
        diagnostic routine d02tzf
        general nonlinear problem solver (thread safe) d02tlf
        interpolation routine d02tyf
        setup routine d02tvf
    finite difference technique with deferred correction, 
        general linear problem d02gbf
        general nonlinear problem, with continuation facility d02raf
        simple nonlinear problem d02gaf
    shooting and matching technique, 
        boundary values to be determined d02haf
        general parameters to be determined d02hbf
        general parameters to be determined, allowing interior matching-point d02agf
        general parameters to be determined, subject to extra algebraic equations d02saf
System of second-order ordinary differential equations, 
    Runge–Kutta–Nystrom method, 
        diagnostic routine d02lyf
        integrator d02laf
        interpolating solutions d02lzf
        setup routine d02lxf

6
Auxiliary Routines Associated with Library Routine Arguments

d02bjwnagf_ode_ivp_rk_zero_simple_dummy_g
See the description of the argument g in d02bjf.
d02bjxnagf_ode_ivp_rk_zero_simple_dummy_output
See the description of the argument output in d02bjf.
d02cjwnagf_ode_ivp_adams_zero_simple_dummy_g
See the description of the argument g in d02cjf.
d02cjxnagf_ode_ivp_adams_zero_simple_dummy_output
See the description of the argument output in d02cjf.
d02ejwnagf_ode_ivp_bdf_zero_simple_dummy_g
See the description of the argument g in d02ejf.
d02ejxnagf_ode_ivp_bdf_zero_simple_dummy_output
See the description of the argument output in d02ejf.
d02ejynagf_ode_ivp_bdf_zero_simple_dummy_pederv
See the description of the argument pederv in d02ejf.
d02gawnagf_ode_bvp_fd_nonlin_gen_dummy_jacobf
See the description of the argument jaceps in d02raf.
d02gaxnagf_ode_bvp_fd_nonlin_gen_dummy_jacobg
See the description of the argument jacgep in d02raf.
d02gaynagf_ode_bvp_fd_nonlin_gen_dummy_jaceps
See the description of the argument jacobg in d02raf.
d02gaznagf_ode_bvp_fd_nonlin_gen_dummy_jacgep
See the description of the argument jacobf in d02raf.
d02hbwnagf_ode_bvp_shoot_genpar_algeq_dummy_prsol
See the description of the argument prsol in d02saf.
d02hbxnagf_ode_bvp_shoot_genpar_algeq_sample_monit
See the description of the argument monit in d02saf.
d02hbynagf_ode_bvp_shoot_genpar_algeq_dummy_constr
See the description of the argument constr in d02saf.
d02hbznagf_ode_bvp_shoot_genpar_algeq_dummy_eqn
See the description of the argument eqn in d02saf.
d02kaynagf_ode_sl2_reg_finite_dummy_monit
See the description of the argument monit in d02kaf, d02kdf and d02kef.
d02nbynagf_ode_ivp_stiff_exp_fulljac_dummy_monit
See the description of the argument monitr in d02nbf, d02ncf, d02ndf, d02ngf, d02nhf and d02njf.
d02nbznagf_ode_ivp_stiff_exp_fulljac_dummy_jac
See the description of the argument jac in d02nbf.
d02ncznagf_ode_ivp_stiff_exp_bandjac_dummy_jac
See the description of the argument jac in d02ncf.
d02ndznagf_ode_ivp_stiff_exp_sparjac_dummy_jac
See the description of the argument jac in d02ndf.
d02neznagf_ode_dae_dassl_gen_dummy_jac
See the description of the argument jac in d02nef.
d02ngznagf_ode_ivp_stiff_imp_fulljac_dummy_jac
See the description of the argument jac in d02ngf.
d02nhznagf_ode_ivp_stiff_imp_bandjac_dummy_jac
See the description of the argument jac in d02nhf.
d02njznagf_ode_ivp_stiff_imp_sparjac_dummy_jac
See the description of the argument jac in d02njf.
d02qfznagf_ode_ivp_adams_roots_dummy_g
See the description of the argument g in d02qff.
d02sasnagf_ode_bvp_shoot_genpar_algeq_dummy_monit
See the description of the argument monit in d02saf.

7
Routines Withdrawn or Scheduled for Withdrawal

The following lists all those routines that have been withdrawn since Mark 19 of the Library or are scheduled for withdrawal at one of the next two marks.
Withdrawn
Routine
Mark of
Withdrawal

Replacement Routine(s)
d02pcf26d02pef and associated D02P routines
d02pdf26d02pff or d02pgf and associated D02P routines
d02pvf26d02pqf
d02pwf26d02prf
d02pxf26d02psf
d02pyf26d02ptf
d02pzf26d02puf
d02tkf27d02tlf

8
References

Ascher U M and Bader G (1987) A new basis implementation for a mixed order boundary value ODE solver SIAM J. Sci. Stat. Comput. 8 483–500
Ascher U M, Christiansen J and Russell R D (1979) A collocation solver for mixed order systems of boundary value problems Math. Comput. 33 659–679
Ascher U M, Mattheij R M M and Russell R D (1988) Numerical Solution of Boundary Value Problems for Ordinary Differential Equations Prentice–Hall
Berzins M, Brankin R W and Gladwell I (1988) Design of the stiff integrators in the NAG Library SIGNUM Newsl. 23 16–23
Brenan K, Campbell S and Petzold L (1996) Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations SIAM, Philadelphia
Gladwell I (1979a) The development of the boundary value codes in the ordinary differential equations chapter of the NAG Library Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer–Verlag
Gladwell I (1979b) Initial value routines in the NAG Library ACM Trans. Math. Software 5 386–400
Gladwell I (1987) The NAG Library boundary value codes Numerical Analysis Report 134 Manchester University
Gladwell I and Sayers D K (ed.) (1980) Computational Techniques for Ordinary Differential Equations Academic Press
Hall G and Watt J M (ed.) (1976) Modern Numerical Methods for Ordinary Differential Equations Clarendon Press, Oxford
Keller H B (1992) Numerical Methods for Two-point Boundary-value Problems Dover, New York
Muite B K (2010) A numerical comparison of Chebyshev methods for solving fourth-order semilinear initial boundary value problems Journal of Computational and Applied Mathematics 234(2) 317–342
Pryce J D (1986) Error estimation for phase-function shooting methods for Sturm–Liouville problems IMA J. Numer. Anal. 6 103–123