# ddriv1.f

          SUBROUTINE DDRIV1 (N, T, Y, F, TOUT, MSTATE, EPS, WORK, LENW,
         8   IERFLG)
    C***BEGIN PROLOGUE  DDRIV1
    C***PURPOSE  The function of DDRIV1 is to solve N (200 or fewer)
    C            ordinary differential equations of the form
    C            dY(I)/dT = F(Y(I),T), given the initial conditions
    C            Y(I) = YI.  DDRIV1 uses double precision arithmetic.
    C***LIBRARY   SLATEC (SDRIVE)
    C***CATEGORY  I1A2, I1A1B
    C***TYPE      DOUBLE PRECISION (SDRIV1-S, DDRIV1-D, CDRIV1-C)
    C***KEYWORDS  DOUBLE PRECISION, GEAR'S METHOD, INITIAL VALUE PROBLEMS,
    C             ODE, ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, STIFF
    C***AUTHOR  Kahaner, D. K., (NIST)
    C             National Institute of Standards and Technology
    C             Gaithersburg, MD  20899
    C           Sutherland, C. D., (LANL)
    C             Mail Stop D466
    C             Los Alamos National Laboratory
    C             Los Alamos, NM  87545
    C***DESCRIPTION
    C
    C   Version 92.1
    C
    C  I.  CHOOSING THE CORRECT ROUTINE  ...................................
    C
    C     SDRIV
    C     DDRIV
    C     CDRIV
    C           These are the generic names for three packages for solving
    C           initial value problems for ordinary differential equations.
    C           SDRIV uses single precision arithmetic.  DDRIV uses double
    C           precision arithmetic.  CDRIV allows complex-valued
    C           differential equations, integrated with respect to a single,
    C           real, independent variable.
    C
    C    As an aid in selecting the proper program, the following is a
    C    discussion of the important options or restrictions associated with
    C    each program:
    C
    C      A. DDRIV1 should be tried first for those routine problems with
    C         no more than 200 differential equations (DDRIV2 and DDRIV3
    C         have no such restriction.)  Internally this routine has two
    C         important technical defaults:
    C           1. Numerical approximation of the Jacobian matrix of the
    C              right hand side is used.
    C           2. The stiff solver option is used.
    C         Most users of DDRIV1 should not have to concern themselves
    C         with these details.
    C
    C      B. DDRIV2 should be considered for those problems for which
    C         DDRIV1 is inadequate.  For example, DDRIV1 may have difficulty
    C         with problems having zero initial conditions and zero
    C         derivatives.  In this case DDRIV2, with an appropriate value
    C         of the parameter EWT, should perform more efficiently.  DDRIV2
    C         provides three important additional options:
    C           1. The nonstiff equation solver (as well as the stiff
    C              solver) is available.
    C           2. The root-finding option is available.
    C           3. The program can dynamically select either the non-stiff
    C              or the stiff methods.
    C         Internally this routine also defaults to the numerical
    C         approximation of the Jacobian matrix of the right hand side.
    C
    C      C. DDRIV3 is the most flexible, and hence the most complex, of
    C         the programs.  Its important additional features include:
    C           1. The ability to exploit band structure in the Jacobian
    C              matrix.
    C           2. The ability to solve some implicit differential
    C              equations, i.e., those having the form:
    C                   A(Y,T)*dY/dT = F(Y,T).
    C           3. The option of integrating in the one step mode.
    C           4. The option of allowing the user to provide a routine
    C              which computes the analytic Jacobian matrix of the right
    C              hand side.
    C           5. The option of allowing the user to provide a routine
    C              which does all the matrix algebra associated with
    C              corrections to the solution components.
    C
    C  II.  PARAMETERS  ....................................................
    C
    C       (REMEMBER--To run DDRIV1 correctly in double precision, ALL
    C       non-integer arguments in the call sequence, including
    C       arrays, MUST be declared double precision.)
    C
    C    The user should use parameter names in the call sequence of DDRIV1
    C    for those quantities whose value may be altered by DDRIV1.  The
    C    parameters in the call sequence are:
    C
    C    N      = (Input) The number of differential equations, N .LE. 200
    C
    C    T      = The independent variable.  On input for the first call, T
    C             is the initial point.  On output, T is the point at which
    C             the solution is given.
    C
    C    Y      = The vector of dependent variables.  Y is used as input on
    C             the first call, to set the initial values.  On output, Y
    C             is the computed solution vector.  This array Y is passed
    C             in the call sequence of the user-provided routine F.  Thus
    C             parameters required by F can be stored in this array in
    C             components N+1 and above.  (Note: Changes by the user to
    C             the first N components of this array will take effect only
    C             after a restart, i.e., after setting MSTATE to +1(-1).)
    C
    C    F      = A subroutine supplied by the user.  The name must be
    C             declared EXTERNAL in the user's calling program.  This
    C             subroutine is of the form:
    C                   SUBROUTINE F (N, T, Y, YDOT)
    C                   DOUBLE PRECISION Y(*), YDOT(*)
    C                     .
    C                     .
    C                   YDOT(1) = ...
    C                     .
    C                     .
    C                   YDOT(N) = ...
    C                   END (Sample)
    C             This computes YDOT = F(Y,T), the right hand side of the
    C             differential equations.  Here Y is a vector of length at
    C             least N.  The actual length of Y is determined by the
    C             user's declaration in the program which calls DDRIV1.
    C             Thus the dimensioning of Y in F, while required by FORTRAN
    C             convention, does not actually allocate any storage.  When
    C             this subroutine is called, the first N components of Y are
    C             intermediate approximations to the solution components.
    C             The user should not alter these values.  Here YDOT is a
    C             vector of length N.  The user should only compute YDOT(I)
    C             for I from 1 to N.  Normally a return from F passes
    C             control back to  DDRIV1.  However, if the user would like
    C             to abort the calculation, i.e., return control to the
    C             program which calls DDRIV1, he should set N to zero.
    C             DDRIV1 will signal this by returning a value of MSTATE
    C             equal to +5(-5).  Altering the value of N in F has no
    C             effect on the value of N in the call sequence of DDRIV1.
    C
    C    TOUT   = (Input) The point at which the solution is desired.
    C
    C    MSTATE = An integer describing the status of integration.  The user
    C             must initialize MSTATE to +1 or -1.  If MSTATE is
    C             positive, the routine will integrate past TOUT and
    C             interpolate the solution.  This is the most efficient
    C             mode.  If MSTATE is negative, the routine will adjust its
    C             internal step to reach TOUT exactly (useful if a
    C             singularity exists beyond TOUT.)  The meaning of the
    C             magnitude of MSTATE:
    C               1  (Input) Means the first call to the routine.  This
    C                  value must be set by the user.  On all subsequent
    C                  calls the value of MSTATE should be tested by the
    C                  user.  Unless DDRIV1 is to be reinitialized, only the
    C                  sign of MSTATE may be changed by the user.  (As a
    C                  convenience to the user who may wish to put out the
    C                  initial conditions, DDRIV1 can be called with
    C                  MSTATE=+1(-1), and TOUT=T.  In this case the program
    C                  will return with MSTATE unchanged, i.e.,
    C                  MSTATE=+1(-1).)
    C               2  (Output) Means a successful integration.  If a normal
    C                  continuation is desired (i.e., a further integration
    C                  in the same direction), simply advance TOUT and call
    C                  again.  All other parameters are automatically set.
    C               3  (Output)(Unsuccessful) Means the integrator has taken
    C                  1000 steps without reaching TOUT.  The user can
    C                  continue the integration by simply calling DDRIV1
    C                  again.
    C               4  (Output)(Unsuccessful) Means too much accuracy has
    C                  been requested.  EPS has been increased to a value
    C                  the program estimates is appropriate.  The user can
    C                  continue the integration by simply calling DDRIV1
    C                  again.
    C               5  (Output)(Unsuccessful) N has been set to zero in
    C                  SUBROUTINE F.
    C               6  (Output)(Successful) For MSTATE negative, T is beyond
    C                  TOUT.  The solution was obtained by interpolation.
    C                  The user can continue the integration by simply
    C                  advancing TOUT and calling DDRIV1 again.
    C               7  (Output)(Unsuccessful) The solution could not be
    C                  obtained.  The value of IERFLG (see description
    C                  below) for a "Recoverable" situation indicates the
    C                  type of difficulty encountered: either an illegal
    C                  value for a parameter or an inability to continue the
    C                  solution.  For this condition the user should take
    C                  corrective action and reset MSTATE to +1(-1) before
    C                  calling DDRIV1 again.  Otherwise the program will
    C                  terminate the run.
    C
    C    EPS    = On input, the requested relative accuracy in all solution
    C             components.  On output, the adjusted relative accuracy if
    C             the input value was too small.  The value of EPS should be
    C             set as large as is reasonable, because the amount of work
    C             done by DDRIV1 increases as EPS decreases.
    C
    C    WORK
    C    LENW   = (Input)
    C             WORK is an array of LENW double precision words used
    C             internally for temporary storage.  The user must allocate
    C             space for this array in the calling program by a statement
    C             such as
    C                       DOUBLE PRECISION WORK(...)
    C             The length of WORK should be at least N*N + 11*N + 300
    C             and LENW should be set to the value used.  The contents of
    C             WORK should not be disturbed between calls to DDRIV1.
    C
    C    IERFLG = An error flag.  The error number associated with a
    C             diagnostic message (see Section IV-A below) is the same as
    C             the corresponding value of IERFLG.  The meaning of IERFLG:
    C               0  The routine completed successfully. (No message is
    C                  issued.)
    C               3  (Warning) The number of steps required to reach TOUT
    C                  exceeds 1000 .
    C               4  (Warning) The value of EPS is too small.
    C              11  (Warning) For MSTATE negative, T is beyond TOUT.
    C                  The solution was obtained by interpolation.
    C              15  (Warning) The integration step size is below the
    C                  roundoff level of T.  (The program issues this
    C                  message as a warning but does not return control to
    C                  the user.)
    C              21  (Recoverable) N is greater than 200 .
    C              22  (Recoverable) N is not positive.
    C              26  (Recoverable) The magnitude of MSTATE is either 0 or
    C                  greater than 7 .
    C              27  (Recoverable) EPS is less than zero.
    C              32  (Recoverable) Insufficient storage has been allocated
    C                  for the WORK array.
    C              41  (Recoverable) The integration step size has gone
    C                  to zero.
    C              42  (Recoverable) The integration step size has been
    C                  reduced about 50 times without advancing the
    C                  solution.  The problem setup may not be correct.
    C             999  (Fatal) The magnitude of MSTATE is 7 .
    C
    C  III.  USAGE  ........................................................
    C
    C                PROGRAM SAMPLE
    C                EXTERNAL F
    C                DOUBLE PRECISION ALFA, EPS, T, TOUT
    C          C                                N is the number of equations
    C                PARAMETER(ALFA = 1.D0, N = 3, LENW = N*N + 11*N + 300)
    C                DOUBLE PRECISION WORK(LENW), Y(N+1)
    C          C                                               Initial point
    C                T = 0.00001D0
    C          C                                      Set initial conditions
    C                Y(1) = 10.D0
    C                Y(2) = 0.D0
    C                Y(3) = 10.D0
    C          C                                              Pass parameter
    C                Y(4) = ALFA
    C                TOUT = T
    C                MSTATE = 1
    C                EPS = .001D0
    C           10   CALL DDRIV1 (N, T, Y, F, TOUT, MSTATE, EPS, WORK, LENW,
    C               8             IERFLG)
    C                IF (MSTATE .GT. 2) STOP
    C                WRITE(*, '(4E12.3)') TOUT, (Y(I), I=1,3)
    C                TOUT = 10.D0*TOUT
    C                IF (TOUT .LT. 50.D0) GO TO 10
    C                END
    C
    C                SUBROUTINE F (N, T, Y, YDOT)
    C                DOUBLE PRECISION ALFA, T, Y(*), YDOT(*)
    C                ALFA = Y(N+1)
    C                YDOT(1) = 1.D0 + ALFA*(Y(2) - Y(1)) - Y(1)*Y(3)
    C                YDOT(2) = ALFA*(Y(1) - Y(2)) - Y(2)*Y(3)
    C                YDOT(3) = 1.D0 - Y(3)*(Y(1) + Y(2))
    C                END
    C
    C  IV.  OTHER COMMUNICATION TO THE USER  ...............................
    C
    C    A. The solver communicates to the user through the parameters
    C       above.  In addition it writes diagnostic messages through the
    C       standard error handling program XERMSG.  A complete description
    C       of XERMSG is given in "Guide to the SLATEC Common Mathematical
    C       Library" by Kirby W. Fong et al..  At installations which do not
    C       have this error handling package the short but serviceable
    C       routine, XERMSG, available with this package, can be used.  That
    C       program uses the file named OUTPUT to transmit messages.
    C
    C    B. The number of evaluations of the right hand side can be found
    C       in the WORK array in the location determined by:
    C            LENW - (N + 50) + 4
    C
    C  V.  REMARKS  ........................................................
    C
    C    For other information, see Section IV of the writeup for DDRIV3.
    C
    C***REFERENCES  C. W. Gear, Numerical Initial Value Problems in
    C                 Ordinary Differential Equations, Prentice-Hall, 1971.
    C***ROUTINES CALLED  DDRIV3, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790601  DATE WRITTEN
    C   900329  Initial submission to SLATEC.
    C***END PROLOGUE  DDRIV1