************************************************************************
*   FUNCTION:  Numerically solve an initial value problem.             *
*   AUTHORS:   Lawrence Shampine, Richard Allen, Steven Pruess  for    *
*              the text  Fundamentals of Numerical Computing           *
*   DATE:      January 3, 1996                                         *
************************************************************************
      SUBROUTINE RKE(F,NEQ,X,Y,H,FIRST,TOL,THRES,FLAG,STEP,YCOEFF,
     &               SCR,NDIM)
C
      INTEGER NEQ,FLAG,NDIM
      LOGICAL FIRST
      DOUBLE PRECISION X,Y(NEQ),H,TOL,THRES(NEQ),STEP,YCOEFF(NDIM),
     &                 SCR(NDIM)
      EXTERNAL F
C
C  RKE integrates a system of NEQ first order ordinary differential
C  equations over one step using a Runge-Kutta method due to R.
C  England. It provides for the evaluation of the solution within
C  the step by interpolation.  RKE is only an interface; the algorithm
C  is implemented in subroutine RKES, the England formulas are
C  evaluated in subroutine ENGLND, and the quintic Hermite interpolant
C  is evaluated in the subroutine YVALUE.  The subroutine YVALUE can
C  be used independently after any successful step by RKE to approx-
C  imate the solution components in the interval [X-STEP,X].
C
C  Input parameters:
C
C   F      = name of the subroutine which defines the system of
C            differential equations.  This subroutine must have
C            the form:
C                 SUBROUTINE F(X,Y,YPRIME)
C                 DOUBLE PRECISION X,Y(*),YPRIME(*)
C                 YPRIME(1) = ...
C                        ...
C                 YPRIME(neq) = ...
C                 RETURN
C                 END
C            The subroutine name F must appear in an EXTERNAL
C            statement in the calling program.
C   NEQ    = number of first order equations in the system.
C   X      = initial value of the independent variable.
C   Y      = vector (with length NEQ) of solution values at X.
C   H      = step size for current step (its sign determines
C            the direction of integration).  On the first call
C            to RKE you must specify H.  After the first call,
C            the code suggests a suitable H for the next step.
C   FIRST  = logical variable indicating first call or
C            subsequent call.  On the first call to RKE set
C            FIRST = .TRUE.; for subsequent calls the code
C            automatically sets FIRST = .FALSE..
C   TOL, THRES = desired tolerances: TOL is a scalar and THRES is a
C            vector with NEQ components.  The value of TOL must
C            be in the interval [ 10*U , 0.01 ] where U is the 
C            unit roundoff for your machine; each component of 
C            THRES must be non-negative.  On the first call, 
C            if some Y(I)=0 then THRES(I) must be positive.  
C            The convergence criterion is roughly
C                 abs(local error in Y(I)) .LE.
C                       TOL*max(THRES(I),abs(Y(I)))
C            for each I.
C   NDIM   = declared dimension of auxiliary storage vector SCR and 
C            the output vector YCOEFF.  NDIM must be at least 6*NEQ.
C
C  Output parameters:
C   X      = value of the independent variable to which the
C            integration advanced.
C   Y      = vector of computed solution values at X.
C   H      = step size suitable for the next step.
C   FLAG   = an integer reporting what the code did:
C              0, if a step was successfully taken;
C              1, if excessive work (more than 500 steps) was
C                 expended;
C              2, if TOL, THRES are too small;
C             -1, if input is invalid.
C   STEP   = output X - input X, the actual length of step
C            taken.
C   YCOEFF = vector (with length 6*NEQ) of coefficient values for the
C            quintic Hermite interpolation used by subroutine YVALUE.
C
C  Auxiliary storage:
C   SCR   = vector of length NDIM (at least 6*NEQ).
C
C  RKE is organized so that subsequent calls to continue the
C  integration involve little, if any, extra effort.  The
C  parameter FLAG must, however, be monitored in order to
C  determine what to do next.  Specifically, if
C      FLAG = 0, the code may be called again to continue the
C                integration another step towards X+H.  You may
C                shorten H if you wish; if its sign is changed,
C                you must restart the integration (FIRST=.TRUE.).
C      FLAG = 1, the code has attempted 500 steps: either TOL/THRES
C                are not appropriate or possibly the solution does 
C                not exist beyond X.  You may alter TOL or THRES.  
C                If you want to continue, set FIRST=.TRUE. and call 
C                RKE again; an additional 500 steps will be allowed.
C      FLAG = 2, the error parameters TOL/THRES are too small.  If
C                a solution component is zero, use a positive value
C                for the corresponding component of THRES.  To
C                continue with increased error parameters, set 
C                FIRST=.TRUE. and call RKE again.
C      FLAG =-1, you cannot continue the solution of this
C                problem; any attempt to do so will result in a
C                program STOP.
C
      IF ((NEQ .LT. 1) .OR. (NDIM .LT. 6*NEQ)) THEN
         FLAG = -1
         RETURN
      ENDIF
      CALL RKES(F,NEQ,X,Y,H,FIRST,TOL,THRES,FLAG,STEP,YCOEFF,NDIM,
     &      SCR(1),SCR(1+NEQ),SCR(1+2*NEQ),SCR(1+3*NEQ),YCOEFF(1),
     &      YCOEFF(1+NEQ),YCOEFF(1+2*NEQ),SCR(1+4*NEQ),YCOEFF(1+3*NEQ),
     &      YCOEFF(1+4*NEQ),YCOEFF(1+5*NEQ),SCR(1+5*NEQ))
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE RKES(F,NEQ,X,Y,H,FIRST,TOL,THRES,FLAG,STEP,YCOEFF,
     &                NDIM,Y4TH,Y5TH,YMID,K0,K1,K2,K3,K4,K5,K6,K7,K8)
      INTEGER NEQ,FLAG,NDIM
      LOGICAL FIRST
      DOUBLE PRECISION X,Y(NEQ),H,TOL,THRES(NEQ),STEP,
     &   YCOEFF(NDIM),Y4TH(NEQ),Y5TH(NEQ),YMID(NEQ),K0(NEQ),
     &   K1(NEQ),K2(NEQ),K3(NEQ),K4(NEQ),K5(NEQ),K6(NEQ),
     &   K7(NEQ),K8(NEQ)
      EXTERNAL F
C
C     Local variables:
C
      INTEGER I,IX,MXSTEP,NASTEP
      LOGICAL FAILED,VALID
      DOUBLE PRECISION ALPHA,BETA,GAMMA,ERR,H12,H180,XMID,TOLAIM,U,WT
      DOUBLE PRECISION ZERO,HNDRTH,TENTH,FIFTH,HALF,SIXTEN,
     &   C1,C2,C4,C6,C8,C10,C12,C14,C15,C24,C32,C64,C92,
     &   C96,C121
      DOUBLE PRECISION C144,C180
      SAVE NASTEP
      DATA ZERO/0.D0/,C1/1.D0/,C2/2.D0/,C4/4.D0/,
     &   C6/6.D0/,C8/8.D0/,C10/10.D0/,C12/12.D0/,C14/14.D0/,
     &   C15/15.D0/,C24/24.D0/,C32/32.D0/,C64/64.D0/,C92/92.D0/,
     &   C96/96.D0/,C121/121.D0/,C144/144.D0/,C180/180.D0/,
     &   HNDRTH/.01D0/,TENTH/.1D0/,FIFTH/.2D0/,HALF/.5D0/,
     &   SIXTEN/.6D0/
      DATA MXSTEP/500/
************************************************************************
*     Machine dependent constant:                                      *
*                                                                      *
*     Set U to the appropriate unit roundoff level for your machine;   *
*     here we assume IEEE double precision.                            *
      DATA U/2.22D-16/
************************************************************************
C
      IF (FIRST) THEN
C
C        First call; test the input data.
C
         VALID = .TRUE.
         IF (NEQ .LT. 1) VALID = .FALSE.
         IF (NDIM .LT. 6*NEQ) VALID = .FALSE.
         IF ((TOL .GT. HNDRTH) .OR.  (TOL .LT. C10*U))
     &      VALID = .FALSE.
         DO 15 I = 1,NEQ
           VALID = VALID .AND. (THRES(I) .GE. ZERO)
   15    CONTINUE
         DO 20 I = 1,NEQ
           IF ((ABS(Y(I)) .EQ. ZERO) .AND. (THRES(I) .EQ. ZERO))
     &               VALID = .FALSE.
   20    CONTINUE
         IF (.NOT. VALID) THEN
            FLAG = -1
            RETURN
         ENDIF
C
C        Initialize.
C
         CALL F(X,Y,K0)
         NASTEP = 0
      ELSE
C
C        Continuation; test input data.
C
         IF (FLAG .EQ. -1) STOP
         VALID = .TRUE.
         IF ((TOL .GT. HNDRTH) .OR. (TOL .LT. C10*U))
     &      VALID = .FALSE.
C
C        Although it is allowed that THRES can be altered between
C        calls, the signs of the components are checked only on the
C        first call.
C
         IF (.NOT. VALID) THEN
            FLAG = -1
            RETURN
         ENDIF
      ENDIF
C
C     Start a new step here.
C
      TOLAIM = SIXTEN*TOL
      FAILED = .FALSE.
C
C     Begin another attempt after a step failure here.
C
   25 IF (NASTEP .GE. MXSTEP) THEN
         FLAG = 1
         RETURN
      ENDIF
C
C     Call ENGLND to take the basic England 4th order step to X+H/2.
C
      CALL ENGLND(F,X,Y,K0,K1,K2,K3,H,YMID,NEQ)
      XMID = X+HALF*H
      CALL F(XMID,YMID,K4)
C
C     Call ENGLND to take basic England 4th order step to X+H.
C
      CALL ENGLND(F,XMID,YMID,K4,K5,K6,K7,H,Y4TH,NEQ)
C
C     Calculate slopes at X+H for the 5th order solution.
C
      H12 = H/C12
      DO 30 I = 1,NEQ
        Y5TH(I) = Y(I)+H12*(-K0(I)-C96*K1(I)+C92*K2(I)
     &        -C121*K3(I)+C144*K4(I)+C6*K5(I)-C12*K6(I))
   30 CONTINUE
      CALL F(X+H,Y5TH,K8)
      H180 = H/C180
      DO 35 I = 1,NEQ
        Y5TH(I) = Y(I)+H180*(C14*K0(I)+C64*K2(I)+C32*K3(I)
     &        -C8*K4(I)+C64*K6(I)+C15*K7(I)-K8(I))
   35 CONTINUE
C
C     Form ERR, the weighted maximum norm of the estimated local error.
C
      ERR = ZERO
      DO 40 I = 1,NEQ
        IF (Y4TH(I) .NE. Y5TH(I)) THEN
           WT = C2*(ABS(Y(I))+ABS(YMID(I)))+ABS(Y4TH(I))+
     &         ABS(Y5TH(I))
           WT = MAX(THRES(I),WT/C6)
           ERR = MAX(ERR,ABS(Y5TH(I)-Y4TH(I))/WT)
        ENDIF
   40 CONTINUE
      NASTEP = NASTEP+1
      IF (ERR .GT. TOL) THEN
C
C        Failed step.
C
         IF (FIRST) THEN
            ALPHA = TENTH
         ELSE
            IF (.NOT. FAILED) THEN
               FAILED = .TRUE.
               ALPHA = (TOLAIM/ERR)**FIFTH
               ALPHA = MAX(TENTH,ALPHA)
            ELSE
               ALPHA = HALF
            ENDIF
         ENDIF
         H = ALPHA*H
         IF (ABS(H) .GE. C24*U*MAX(ABS(X),ABS(X+H))) GOTO 25
         FLAG = 2
      ELSE
C
C        Successful step.
C
         X = X+H
         CALL F(X,Y5TH,K8)
         STEP = H
         DO 45 I = 1,NEQ
            ALPHA = (K0(I)-K8(I))*STEP
            BETA = Y5TH(I)-C2*YMID(I)+Y(I)
            GAMMA = Y5TH(I)-Y(I)
            IX = 6*(I-1)
            YCOEFF(1+IX) = YMID(I)
            YCOEFF(2+IX) = K4(I)*STEP
            YCOEFF(3+IX) = C4*BETA+HALF*ALPHA
            YCOEFF(4+IX) = C10*GAMMA-STEP*(K0(I)+C8*K4(I)+K8(I))
            YCOEFF(5+IX) = -C8*BETA-C2*ALPHA
            YCOEFF(6+IX) = -C24*GAMMA+C4*STEP*
     &                      (K0(I)+C4*K4(I)+K8(I))
            Y(I) = Y5TH(I)
            K0(I) = K8(I)
   45    CONTINUE
         FIRST = .FALSE.
         BETA = (ERR/TOLAIM)**FIFTH
         ALPHA = C1/MAX(TENTH,BETA)
         IF (FAILED) ALPHA = MIN(C1,ALPHA)
         H = ALPHA*H
         FLAG = 0
      ENDIF
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE ENGLND(F,X,Y,FIN,FOUT1,FOUT2,FOUT3,H,YNEW,NEQ)
      INTEGER NEQ,I
      DOUBLE PRECISION X,Y(NEQ),FIN(NEQ),FOUT1(NEQ),FOUT2(NEQ),
     &   FOUT3(NEQ),H,YNEW(NEQ)
      DOUBLE PRECISION H2,H4,H8,H12,EIGHTH,QUART,HALF,TWO,
     &   FOUR,TWELVE
      EXTERNAL F
      DATA EIGHTH/0.125D0/,QUART/0.25D0/,HALF/0.5D0/,TWO/2.D0/
      DATA FOUR/4.D0/,TWELVE/12.D0/
C
C     Use the basic England formulas to compute YNEW, the 4th
C     order solution at X+H/2.
C
      H4 = QUART*H
      DO 10 I = 1,NEQ
        YNEW(I) = Y(I)+H4*FIN(I)
  10  CONTINUE
      CALL F(X+H4,YNEW,FOUT1)
C
      H8 = EIGHTH*H
      DO 20 I = 1,NEQ
        YNEW(I) = Y(I)+H8*(FIN(I)+FOUT1(I))
  20  CONTINUE
      CALL F(X+H4,YNEW,FOUT2)
C
      H2 = HALF*H
      DO 30 I = 1,NEQ
        YNEW(I) = Y(I)+H2*(-FOUT1(I)+TWO*FOUT2(I))
  30  CONTINUE
      CALL F(X+H2,YNEW,FOUT3)
C
      H12 = H/TWELVE
      DO 40 I = 1,NEQ
        YNEW(I) = Y(I)+H12*(FIN(I)+FOUR*FOUT2(I)+FOUT3(I))
  40  CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE YVALUE(NEQ,X,STEP,YCOEFF,POINT,U)
      INTEGER NEQ
      DOUBLE PRECISION X,STEP,YCOEFF(*),POINT,U(NEQ)
C
C     Evaluate the quintic Hermite interpolants based on output from
C     subprogram RKE.
C
C  Input parameters:
C    NEQ    =     number of components in vector Y.
C    X      =     output from RKE (final value for independent
C                 variable).
C    STEP   =     output from RKE (total step taken).
C    YCOEFF = output from RKE (coefficients of quintics).
C    POINT  =     point at which evaluation is to be made.
C  Output parameters:
C    U      =     vector of solution components at POINT.
C
C  Local variables
      INTEGER I,IX
      DOUBLE PRECISION Z,HALF
      DATA HALF/0.5D0/
C
C     Transform from [X-STEP,X] to [-0.5,0.5], then evaluate.
C
      Z = (POINT-X)/STEP+HALF
      DO 10 I = 1,NEQ
        IX = 6*(I-1)
        U(I) = YCOEFF(1+IX)+Z*(YCOEFF(2+IX)+Z*(YCOEFF(3+IX)+Z*
     &        (YCOEFF(4+IX)+Z*(YCOEFF(5+IX)+Z*YCOEFF(6+IX)))))
   10 CONTINUE
      RETURN
      END
C=======================================================================
