# dsteps.f

          SUBROUTINE DSTEPS (DF, NEQN, Y, X, H, EPS, WT, START, HOLD, K,
         +   KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G,
         +   PHASE1, NS, NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV,
         +   KGI, GI, RPAR, IPAR)
    C***BEGIN PROLOGUE  DSTEPS
    C***PURPOSE  Integrate a system of first order ordinary differential
    C            equations one step.
    C***LIBRARY   SLATEC (DEPAC)
    C***CATEGORY  I1A1B
    C***TYPE      DOUBLE PRECISION (STEPS-S, DSTEPS-D)
    C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
    C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
    C***AUTHOR  Shampine, L. F., (SNLA)
    C           Gordon, M. K., (SNLA)
    C             MODIFIED BY H.A. WATTS
    C***DESCRIPTION
    C
    C   Written by L. F. Shampine and M. K. Gordon
    C
    C   Abstract
    C
    C   Subroutine  DSTEPS  is normally used indirectly through subroutine
    C   DDEABM .  Because  DDEABM  suffices for most problems and is much
    C   easier to use, using it should be considered before using  DSTEPS
    C   alone.
    C
    C   Subroutine DSTEPS integrates a system of  NEQN  first order ordinary
    C   differential equations one step, normally from X to X+H, using a
    C   modified divided difference form of the Adams Pece formulas.  Local
    C   extrapolation is used to improve absolute stability and accuracy.
    C   The code adjusts its order and step size to control the local error
    C   per unit step in a generalized sense.  Special devices are included
    C   to control roundoff error and to detect when the user is requesting
    C   too much accuracy.
    C
    C   This code is completely explained and documented in the text,
    C   Computer Solution of Ordinary Differential Equations, The Initial
    C   Value Problem  by L. F. Shampine and M. K. Gordon.
    C   Further details on use of this code are available in "Solving
    C   Ordinary Differential Equations with ODE, STEP, and INTRP",
    C   by L. F. Shampine and M. K. Gordon, SLA-73-1060.
    C
    C
    C   The parameters represent --
    C      DF -- subroutine to evaluate derivatives
    C      NEQN -- number of equations to be integrated
    C      Y(*) -- solution vector at X
    C      X -- independent variable
    C      H -- appropriate step size for next step.  Normally determined by
    C           code
    C      EPS -- local error tolerance
    C      WT(*) -- vector of weights for error criterion
    C      START -- logical variable set .TRUE. for first step,  .FALSE.
    C           otherwise
    C      HOLD -- step size used for last successful step
    C      K -- appropriate order for next step (determined by code)
    C      KOLD -- order used for last successful step
    C      CRASH -- logical variable set .TRUE. when no step can be taken,
    C           .FALSE. otherwise.
    C      YP(*) -- derivative of solution vector at  X  after successful
    C           step
    C      KSTEPS -- counter on attempted steps
    C      TWOU -- 2.*U where U is machine unit roundoff quantity
    C      FOURU -- 4.*U where U is machine unit roundoff quantity
    C      RPAR,IPAR -- parameter arrays which you may choose to use
    C            for communication between your program and subroutine F.
    C            They are not altered or used by DSTEPS.
    C   The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G,
    C   W,P,IV and GI are required for the interpolation subroutine SINTRP.
    C   The remaining variables and arrays are included in the call list
    C   only to eliminate local retention of variables between calls.
    C
    C   Input to DSTEPS
    C
    C      First call --
    C
    C   The user must provide storage in his calling program for all arrays
    C   in the call list, namely
    C
    C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
    C    1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
    C    2  RPAR(*),IPAR(*)
    C
    C    **Note**
    C
    C   The user must also declare  START ,  CRASH ,  PHASE1  and  NORND
    C   logical variables and  DF  an EXTERNAL subroutine, supply the
    C   subroutine  DF(X,Y,YP)  to evaluate
    C      DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN))
    C   and initialize only the following parameters.
    C      NEQN -- number of equations to be integrated
    C      Y(*) -- vector of initial values of dependent variables
    C      X -- initial value of the independent variable
    C      H -- nominal step size indicating direction of integration
    C           and maximum size of step.  Must be variable
    C      EPS -- local error tolerance per step.  Must be variable
    C      WT(*) -- vector of non-zero weights for error criterion
    C      START -- .TRUE.
    C      YP(*) -- vector of initial derivative values
    C      KSTEPS -- set KSTEPS to zero
    C      TWOU -- 2.*U where U is machine unit roundoff quantity
    C      FOURU -- 4.*U where U is machine unit roundoff quantity
    C   Define U to be the machine unit roundoff quantity by calling
    C   the function routine  D1MACH,  U = D1MACH(4), or by
    C   computing U so that U is the smallest positive number such
    C   that 1.0+U .GT. 1.0.
    C
    C   DSTEPS  requires that the L2 norm of the vector with components
    C   LOCAL ERROR(L)/WT(L)  be less than  EPS  for a successful step.  The
    C   array  WT  allows the user to specify an error test appropriate
    C   for his problem.  For example,
    C      WT(L) = 1.0  specifies absolute error,
    C            = ABS(Y(L))  error relative to the most recent value of the
    C                 L-th component of the solution,
    C            = ABS(YP(L))  error relative to the most recent value of
    C                 the L-th component of the derivative,
    C            = MAX(WT(L),ABS(Y(L)))  error relative to the largest
    C                 magnitude of L-th component obtained so far,
    C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  specifies a mixed
    C                 relative-absolute test where  RELERR  is relative
    C                 error,  ABSERR  is absolute error and  EPS =
    C                 MAX(RELERR,ABSERR) .
    C
    C      Subsequent calls --
    C
    C   Subroutine  DSTEPS  is designed so that all information needed to
    C   continue the integration, including the step size  H  and the order
    C   K , is returned with each step.  With the exception of the step
    C   size, the error tolerance, and the weights, none of the parameters
    C   should be altered.  The array  WT  must be updated after each step
    C   to maintain relative error tests like those above.  Normally the
    C   integration is continued just beyond the desired endpoint and the
    C   solution interpolated there with subroutine  SINTRP .  If it is
    C   impossible to integrate beyond the endpoint, the step size may be
    C   reduced to hit the endpoint since the code will not take a step
    C   larger than the  H  input.  Changing the direction of integration,
    C   i.e., the sign of  H , requires the user set  START = .TRUE. before
    C   calling  DSTEPS  again.  This is the only situation in which  START
    C   should be altered.
    C
    C   Output from DSTEPS
    C
    C      Successful Step --
    C
    C   The subroutine returns after each successful step with  START  and
    C   CRASH  set .FALSE. .  X  represents the independent variable
    C   advanced one step of length  HOLD  from its value on input and  Y
    C   the solution vector at the new value of  X .  All other parameters
    C   represent information corresponding to the new  X  needed to
    C   continue the integration.
    C
    C      Unsuccessful Step --
    C
    C   When the error tolerance is too small for the machine precision,
    C   the subroutine returns without taking a step and  CRASH = .TRUE. .
    C   An appropriate step size and error tolerance for continuing are
    C   estimated and all other information is restored as upon input
    C   before returning.  To continue with the larger tolerance, the user
    C   just calls the code again.  A restart is neither required nor
    C   desirable.
    C
    C***REFERENCES  L. F. Shampine and M. K. Gordon, Solving ordinary
    C                 differential equations with ODE, STEP, and INTRP,
    C                 Report SLA-73-1060, Sandia Laboratories, 1973.
    C***ROUTINES CALLED  D1MACH, DHSTRT
    C***REVISION HISTORY  (YYMMDD)
    C   740101  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   890831  Modified array declarations.  (WRB)
    C   890831  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  DSTEPS