      SUBROUTINE SETUP(NEQ,TSTART,YSTART,TEND,TOL,THRES,METHOD,TASK,
     &                 ERRASS,HSTART,WORK,LENWRK,MESAGE)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code SETUP and how it is used in
C  conjunction with UT or CT to solve initial value problems, you should study
C  the document file rksuite.doc carefully before attempting to use the code.
C  The following "Brief Reminder" is intended only to remind you of the
C  meaning, type, and size requirements of the arguments.
C
C  The environmental parameters OUTCH, MCHEPS, and DWARF are used in the
C  following description.  To find out their values
C
C       CALL ENVIRN(OUTCH,MCHEPS,DWARF)
C
C  INPUT VARIABLES
C
C     NEQ       - INTEGER
C                 The number of differential equations in the system.
C                 Constraint: NEQ >= 1
C     TSTART    - DOUBLE PRECISION
C                 The initial value of the independent variable.
C     YSTART(*) - DOUBLE PRECISION array of length NEQ
C                 The vector of initial values of the solution components.
C     TEND      - DOUBLE PRECISION
C                 The integration proceeds from TSTART in the direction of
C                 TEND. You cannot go past TEND.
C                 Constraint: TEND must be clearly distinguishable from TSTART
C                 in the precision available.
C     TOL       - DOUBLE PRECISION
C                 The relative error tolerance.
C                 Constraint: 0.01D0 >= TOL >= 10*MCHEPS
C     THRES(*)  - DOUBLE PRECISION array of length NEQ
C                 THRES(L) is the threshold for the Ith solution component.
C                 Constraint: THRES(L) >= SQRT(DWARF)
C     METHOD    - INTEGER
C                 Specifies which Runge-Kutta pair is to be used.
C                  = 1 - use the (2,3) pair
C                  = 2 - use the (4,5) pair
C                  = 3 - use the (7,8) pair
C     TASK      - CHARACTER*(*)
C                 Only the first character of TASK is significant.
C                 TASK(1:1) = `U' or `u' - UT is to be used
C                           = `C' or `c' - CT is to be used
C                 Constraint: TASK(1:1) = `U'or `u' or`C' or `c'
C     ERRASS    - LOGICAL
C                 = .FALSE. - do not attempt to assess the true error.
C                 = .TRUE.  - assess the true error. Costs roughly twice
C                             as much as the integration with METHODs 2 and
C                             3, and three times with METHOD = 1.
C     HSTART    - DOUBLE PRECISION
C                 0.0D0     - select automatically the first step size.
C                 non-zero  - try HSTART for the first step.
C
C  WORKSPACE
C
C     WORK(*) - DOUBLE PRECISION array of length LENWRK
C               Do not alter the contents of this array after calling SETUP.
C
C  INPUT VARIABLES
C
C     LENWRK  - INTEGER
C               Length of WORK(*): How big LENWRK must be depends
C               on the task and how it is to be solved.
C
C               LENWRK = 32*NEQ is sufficient for all cases. 
C
C               If storage is a problem, the least storage possible 
C               in the various cases is:
C
C                 If TASK = `U' or `u', then
C                   if ERRASS = .FALSE. and
C                     METHOD = 1, LENWRK must be at least 10*NEQ
C                            = 2                          20*NEQ
C                            = 3                          16*NEQ
C                   if ERRASS = .TRUE. and
C                     METHOD = 1, LENWRK must be at least 15*NEQ
C                            = 2                          32*NEQ
C                            = 3                          21*NEQ
C
C                 If TASK = `C' or `c', then
C                   if ERRASS = .FALSE. and
C                     METHOD = 1, LENWRK must be at least 10*NEQ
C                            = 2                          14*NEQ
C                            = 3                          16*NEQ
C                   if ERRASS = .TRUE. and
C                     METHOD = 1, LENWRK must be at least 15*NEQ
C                            = 2                          26*NEQ
C                            = 3                          21*NEQ
C
C                 Warning:  To exploit the interpolation capability
C                 of METHODs 1 and 2, you have to call INTRP.  This
C                 subroutine requires working storage in addition to
C                 that specified here.
C
C     MESAGE    - LOGICAL
C                 Specifies whether you want informative messages written to
C                 the standard output channel OUTCH.
C                 = .TRUE.   - provide messages
C                 = .FALSE.  - do not provide messages
C
C  In the event of a "catastrophic" failure to call SETUP correctly, the
C  nature of the catastrophe is reported on the standard output channel,
C  regardless of the value of MESAGE.  Unless special provision was made
C  in advance (see rksuite.doc), the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HSTART, TEND, TOL, TSTART
      INTEGER           LENWRK, METHOD, NEQ
      LOGICAL           ERRASS, MESAGE
      CHARACTER*(*)     TASK
C     .. Array Arguments ..
      DOUBLE PRECISION  THRES(*), WORK(*), YSTART(*)
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block for General Workspace Pointers ..
      INTEGER           PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      COMMON /RKCOM3/   PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      SAVE   /RKCOM3/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='SETUP')
      INTEGER           MINUS1
      LOGICAL           TELL
      PARAMETER         (MINUS1=-1,TELL=.FALSE.)
      DOUBLE PRECISION  ONE, ZERO, PT01
      PARAMETER         (ONE=1.0D+0,ZERO=0.0D+0,PT01=0.01D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  HMIN
      INTEGER           FLAG, FREEPR, IER, L, LINTPL, LREQ, NREC, VECSTG
      LOGICAL           LEGALT, REQSTG
      CHARACTER         TASK1
C     .. External Subroutines ..
      EXTERNAL          CONST, MCONST, RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SIGN
C     .. Executable Statements ..
C
C  Clear previous flag values of subprograms in the suite.
C
      IER = MINUS1
      CALL RKSIT(TELL,SRNAME,IER)
C
      IER = 1
      NREC = 0
C
C  Fetch output channel and machine constants; initialise common
C  block /RKCOM7/
C
      CALL MCONST(METHOD)
C
C  Check for valid input of trivial arguments
      TASK1 = TASK(1:1)
      LEGALT = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u' .OR.
     &         TASK1 .EQ. 'C' .OR. TASK1 .EQ. 'c'
      IF (.NOT.LEGALT) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A,A,A/A)')
     &' ** You have set the first character of ',
     &' ** TASK to be ''',TASK1,'''. It must be one of ',
     &' ** ''U'',''u'',''C'' or ''c''.'
      ELSE IF (NEQ.LT.1) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A,I6,A)')
     &' ** You have set NEQ = ',NEQ,' which is less than 1.'
      ELSE IF (METHOD.LT.1 .OR. METHOD.GT.3) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A,I6,A)')
     &' ** You have set METHOD = ',METHOD,' which is not 1, 2, or 3.'
      ELSE IF (TSTART.EQ.TEND) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A,D13.5,A)')
     &' ** You have set TSTART = TEND = ',TSTART,'.'
      ELSE IF ((TOL.GT.PT01) .OR. (TOL.LT.RNDOFF)) THEN
         IER = 911
         NREC = 2
         WRITE (REC,'(A,D13.5,A/A,D13.5,A)')
     &' ** You have set TOL = ',TOL,' which is not permitted. The',
     &' ** range of permitted values is (',RNDOFF,',0.01D0).'
      ELSE
         L = 1
   20    CONTINUE
         IF (THRES(L).LT.TINY) THEN
            IER = 911
            NREC = 2
            WRITE (REC,'(A,I6,A,D13.5,A/A,D13.5,A)')
     &' ** You have set THRES(',L,') to be ',THRES(L),' which is ',
     &' ** less than the permitted minimum,',TINY,'.'
         END IF
         L = L + 1
         IF (IER.NE.911 .AND. L.LE.NEQ) GO TO 20
      END IF
C
C  Return if error detected
C
      IF (IER.NE.1) GO TO 80
C
C  Set formula definitions and characteristics by means of arguments
C  in the call list and COMMON blocks /RKCOM4/ and /RKCOM5/
C
      CALL CONST(METHOD,VECSTG,REQSTG,LINTPL)
C
C  Set options in /RKCOM8/
      UTASK = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u'
      MSG = MESAGE
C
C  Initialise problem status in /RKCOM1/ and /RKCOM2/
      NEQN = NEQ
      TSTRT = TSTART
      TND = TEND
      T = TSTART
      TOLD = TSTART
      DIR = SIGN(ONE,TEND-TSTART)
C
C  In CT the first step taken will have magnitude H.  If HSTRT = ABS(HSTART)
C  is not equal to zero, H = HSTRT.  If HSTRT is equal to zero, the code is
C  to find an on-scale initial step size H.  To start this process, H is set
C  here to an upper bound on the first step size that reflects the scale of
C  the independent variable.  UT has some additional information, namely the
C  first output point, that is used to refine this bound in UT when UTASK
C  is .TRUE..  If HSTRT is not zero, but it is either too big or too small,
C  the input HSTART is ignored and HSTRT is set to zero to activate the
C  automatic determination of an on-scale initial step size.
C
      HSTRT = ABS(HSTART)
      HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND)))
      IF (HSTRT.GT.ABS(TEND-TSTART) .OR. HSTRT.LT.HMIN) HSTRT = ZERO
      IF (HSTRT.EQ.ZERO) THEN
         H = MAX(ABS(TEND-TSTART)/RS3,HMIN)
      ELSE
         H = HSTRT
      END IF
      HOLD = ZERO
      TOLR = TOL
      NFCN = 0
      SVNFCN = 0
      OKSTP = 0
      FLSTP = 0
      FIRST = .TRUE.
      LAST = .FALSE.
C
C  WORK(*) is partioned into a number of arrays using pointers. These
C  pointers are set in /RKCOM3/.
      PRTHRS = 1
C                           the threshold values
      PRERST = PRTHRS + NEQ
C                           the error estimates
      PRWT = PRERST + NEQ
C                           the weights used in the local error test
      PRYOLD = PRWT + NEQ
C                           the previous value of the solution
      PRSCR = PRYOLD + NEQ
C                           scratch array used for the higher order
C                           approximate solution and for the previous 
C                           value of the derivative of the solution
      PRY = PRSCR + NEQ
C                           the dependent variables
      PRYP = PRY + NEQ
C                           the derivatives
      PRSTGS = PRYP + NEQ
C                           intermediate stages held in an internal
C                           array STAGES(NEQ,VECSTG)
C
      FREEPR = PRSTGS + VECSTG*NEQ
C
C  Allocate storage for interpolation if the TASK = `U' or `u' was
C  specified. INTP and LINTPL returned by CONST indicate whether there
C  is an interpolation scheme associated with the pair and how much
C  storage is required.
C
      PRINTP = 1
      LNINTP = 1
      IF (UTASK) THEN
         IF (INTP) THEN
            LNINTP = LINTPL*NEQ
            IF (REQSTG) THEN
               PRINTP = FREEPR
               FREEPR = PRINTP + LNINTP
            ELSE
               PRINTP = PRSTGS
               FREEPR = MAX(PRINTP+VECSTG*NEQ,PRINTP+LNINTP)
            END IF
         END IF
      END IF
C
C  Initialise state and allocate storage for global error assessment
C  using /RKCOM6/
      GNFCN = 0
      MAXERR = ZERO
      LOCMAX = TSTART
      ERASON = ERRASS
      ERASFL = .FALSE.
      IF (ERRASS) THEN
C
C  Storage is required for the stages of a secondary integration. The
C  stages of the primary intergration can only be overwritten in the
C  cases where there is no interpolant or the interpolant does not
C  require information about the stages (e.g. METHOD 3 and METHOD 1,
C  respectively).
         IF (.NOT.REQSTG) THEN
            PRZSTG = PRSTGS
         ELSE
            PRZSTG = FREEPR
            FREEPR = PRZSTG + VECSTG*NEQ
         END IF
         PRZY = FREEPR
         PRZYP = PRZY + NEQ
         PRZERS = PRZYP + NEQ
         PRZERR = PRZERS + NEQ
         PRZYNU = PRZERR + NEQ
         FREEPR = PRZYNU + NEQ
      ELSE
         PRZSTG = 1
         PRZY = 1
         PRZYP = 1
         PRZERS = 1
         PRZERR = 1
         PRZYNU = 1
      END IF
C
      LREQ = FREEPR - 1
C
C  Check for enough workspace and suitable range of integration
C
      IF (LENWRK.LT.LREQ) THEN
         IER = 911
         NREC = 2
         WRITE (REC,'(A/A,I6,A,I6,A)')
     &' ** You have not supplied enough workspace. You gave LENWRK ',
     &' ** as', LENWRK, ', but it must be at least ',LREQ,'.'
      ELSE
         HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND)))
         IF (ABS(TEND-TSTART).LT.HMIN) THEN
            IER = 911
            NREC = 4
            WRITE (REC,'(A/A/A,D13.5/A,D13.5,A)')
     &' ** You have set values for TEND and TSTART that are not ',
     &' ** clearly distinguishable for the method and the precision ',
     &' ** of the computer being used. ABS(TEND-TSTART) is ',
     &ABS(TEND-TSTART),
     &' ** but should be at least ',HMIN,'.'
         END IF
      END IF
C
C  Return if error detected
C
      IF (IER.NE.1) GO TO 80
C
C  Initialize elements of the workspace
      DO 40 L = 1, NEQ
         WORK(PRTHRS-1+L) = THRES(L)
         WORK(PRY-1+L) = YSTART(L)
   40 CONTINUE
C
C  Initialize the global error to zero when ERRASS = .TRUE.
      IF (ERRASS) THEN
         DO 60 L = 1, NEQ
            WORK(PRZERR-1+L) = ZERO
   60    CONTINUE
      END IF
C
   80 CONTINUE
C
      CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
      RETURN
      END
      SUBROUTINE UT(F,TWANT,TGOT,YGOT,YPGOT,YMAX,WORK,UFLAG)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code UT and how it is used in
C  conjunction with SETUP to solve initial value problems, you should study
C  the document file rksuite.doc carefully before proceeding further.  The
C  following "Brief Reminder" is intended only to remind you of the meaning,
C  type, and size requirements of the arguments.
C
C  NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM:
C
C     F         - name of the subroutine for evaluating the differential
C                 equations.
C
C  The subroutine F must have the form
C
C  SUBROUTINE F(T,Y,YP)
C  DOUBLE PRECISION T,Y(*),YP(*)
C     Given input values of the independent variable T and the solution
C     components Y(*), for each L = 1,2,...,NEQ evaluate the differential
C     equation for the derivative of the Ith solution component and place the
C     value in YP(L).  Do not alter the input values of T and Y(*).
C  RETURN
C  END
C
C  INPUT VARIABLE
C
C     TWANT     - DOUBLE PRECISION
C                 The next value of the independent variable where a
C                 solution is desired.
C
C                 Constraints: TWANT must lie between the previous value
C                 of TGOT (TSTART on the first call) and TEND. TWANT can be
C                 equal to TEND, but it must be clearly distinguishable from
C                 the previous value of TGOT (TSTART on the first call) in 
C                 the precision available.
C
C  OUTPUT VARIABLES
C
C     TGOT      - DOUBLE PRECISION
C                 A solution has been computed at this value of the
C                 independent variable.
C     YGOT(*)   - DOUBLE PRECISION array of length NEQ
C                 Approximation to the true solution at TGOT. Do not alter
C                 the contents of this array
C     YPGOT(*)  - DOUBLE PRECISION array of length NEQ
C                 Approximation to the first derivative of the true
C                 solution at TGOT.
C     YMAX(*)   - DOUBLE PRECISION array of length NEQ
C                 YMAX(L) is the largest magnitude of YGOT(L) computed at any
C                 time in the integration from TSTART to TGOT. Do not alter
C                 the contents of this array.
C
C  WORKSPACE
C
C     WORK(*)   - DOUBLE PRECISION array as used in SETUP
C                 Do not alter the contents of this array.
C
C  OUTPUT VARIABLE
C
C     UFLAG     - INTEGER
C
C                       SUCCESS.  TGOT = TWANT.
C                 = 1 - Complete success.
C
C                       "SOFT" FAILURES
C                 = 2 - Warning:  You are using METHOD = 3 inefficiently
C                       by computing answers at many values of TWANT.  If
C                       you really need answers at so many specific points,
C                       it would be more efficient to compute them with
C                       METHOD = 2.  To do this you would need to restart
C                       from TGOT, YGOT(*) by a call to SETUP.  If you wish
C                       to continue as you are, you may.
C                 = 3 - Warning:  A considerable amount of work has been
C                       expended.  If you wish to continue on to TWANT, just
C                       call UT again.
C                 = 4 - Warning:  It appears that this problem is "stiff".
C                       You really should change to another code that is
C                       intended for such problems, but if you insist, you can
C                       continue with UT by calling it again.
C
C                       "HARD" FAILURES
C                 = 5 - You are asking for too much accuracy. You cannot
C                       continue integrating this problem.
C                 = 6 - The global error assessment may not be reliable beyond
C                       the current point in the integration.  You cannot
C                       continue integrating this problem.
C
C                       "CATASTROPHIC" FAILURES
C                 = 911 - The nature of the catastrophe is reported on
C                         the standard output channel. Unless special
C                         provision was made in advance (see rksuite.doc),
C                         the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TGOT, TWANT
      INTEGER           UFLAG
C     .. Array Arguments ..
      DOUBLE PRECISION  WORK(*), YGOT(*), YMAX(*), YPGOT(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block for General Workspace Pointers ..
      INTEGER           PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      COMMON /RKCOM3/   PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      SAVE   /RKCOM3/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='UT')
      LOGICAL           ASK, TELL
      PARAMETER         (ASK=.TRUE.,TELL=.FALSE.)
      INTEGER           MINUS1, MINUS2
      PARAMETER         (MINUS1=-1,MINUS2=-2)
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  HMIN, TLAST, TNOW, UTEND
      INTEGER           CFLAG, IER, L, NREC, STATE
      LOGICAL           BADERR, GOBACK
C     .. External Subroutines ..
      EXTERNAL          CHKFL, CT, INTRP, RESET, RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MIN
C     .. Save statement ..
      SAVE              UTEND, TLAST
C     .. Executable Statements ..
      IER = 1
      NREC = 0
      GOBACK = .FALSE.
      BADERR = .FALSE.
C
C  Is it permissible to call UT?
C
      CALL RKSIT(ASK,'SETUP',STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 100
      END IF
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A)')
     &' ** You have not called SETUP, so you cannot use UT.'
         GO TO 100
      END IF
      IF (.NOT.UTASK) THEN
         IER = 911
         NREC = 2
         WRITE (REC,'(A/A)')
     &' ** You have called UT after you specified in SETUP that ',
     &' ** you were going to use CT. This is not permitted.'
         GO TO 100
      END IF
      CALL RKSIT(ASK,SRNAME,STATE)
      IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A/A)')
     &' ** This routine has already returned with a hard failure.',
     &' ** You must call SETUP to start another problem.'
         GO TO 100
      END IF
      STATE = MINUS2
      CALL RKSIT(TELL,SRNAME,STATE)
C
      IF (FIRST) THEN
C
C  First call.
C
C  A value of TND is specified in SETUP. When INTP = .FALSE., as with
C  METHD = 3, output is obtained at the specified TWANT by resetting TND
C  to TWANT.  At this point, before the integration gets started, this can
C  be done with a simple assignment.  Later it is done with a call to RESET.
C  The original TND is SAVEd as a local variable UTEND.
C
         UTEND = TND
         IF (.NOT.INTP) TND = TWANT
C
C  The last TGOT returned is SAVEd in the variable TLAST.  T (a variable
C  passed through the common block RKCOM2) records how far the integration
C  has advanced towards the specified TND.  When output is obtained by
C  interpolation, the integration goes past the TGOT returned (T is closer
C  to the specified TND than TGOT).  Initialize these variables and YMAX(*).
         TLAST = TSTRT
         TGOT = TSTRT
         DO 20 L = 1, NEQN
            YMAX(L) = ABS(WORK(PRY-1+L))
   20    CONTINUE
C
C  If the code is to find an on-scale initial step size H, a bound was placed
C  on H in SETUP.  Here the first output point is used to refine this bound.
         IF (HSTRT.EQ.ZERO) THEN
            H = MIN(ABS(H),ABS(TWANT-TSTRT))
            HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND)))
            H = MAX(H,HMIN)
         END IF
C
      ELSE
C
C  Subsequent call.
C
         IF (TLAST.EQ.UTEND) THEN
            IER = 911
            NREC = 3
            WRITE (REC,'(A/A/A)')
     &' ** You have called UT after reaching TEND. (Your last    ',
     &' ** call to UT resulted in TGOT = TEND.)  To start a new ',
     &' ** problem, you will need to call SETUP.'
            GO TO 100
         END IF
C
      END IF
C
C  Check for valid TWANT.
C
      IF (DIR*(TWANT-TLAST).LE.ZERO) THEN
         IER = 911
         NREC = 4
         WRITE (REC,'(A/A/A/A)')
     &' ** You have made a call to UT with a TWANT that does   ',
     &' ** not lie between the previous value of TGOT (TSTART  ',
     &' ** on the first call) and TEND. This is not permitted. ',
     &' ** Check your program carefully.'
         GO TO 100
      END IF
      IF (DIR*(TWANT-UTEND).GT.ZERO) THEN
         HMIN = MAX(TINY,TOOSML*MAX(ABS(TWANT),ABS(UTEND)))
         IF (ABS(TWANT-UTEND).LT.HMIN) THEN
            IER = 911
            NREC = 5
            WRITE (REC,'(A/A/A/A)')
     &' ** You have made a call to UT with a TWANT that does      ',
     &' ** not lie between the previous value of TGOT (TSTART on  ',
     &' ** the first call) and TEND. This is not permitted. TWANT ',
     &' ** is very close to TEND, so you may have meant to set    ',
     &' ** it to be TEND exactly.  Check your program carefully.  '
         ELSE
            IER = 911
            NREC = 4
            WRITE (REC,'(A/A/A/A)')
     &' ** You have made a call to UT with a TWANT that does   ',
     &' ** not lie between the previous value of TGOT (TSTART  ',
     &' ** on the first call) and TEND. This is not permitted. ',
     &' ** Check your program carefully.'
         END IF
         GO TO 100
      END IF
      IF (.NOT.INTP) THEN
         HMIN = MAX(TINY,TOOSML*MAX(ABS(TLAST),ABS(TWANT)))
         IF (ABS(TWANT-TLAST).LT.HMIN) THEN
            IER = 911
            NREC = 4
            WRITE (REC,'(A/A/A/A,D13.5,A)')
     &' ** You have made a call to UT with a TWANT that is not ',
     &' ** sufficiently different from the last value of TGOT  ',
     &' ** (TSTART on the first call).  When using METHOD = 3, ',
     &' ** it must differ by at least ',HMIN,'.'
            GO TO 100
         END IF
C
C  We have a valid TWANT. There is no interpolation with this METHD and
C  therefore we step to TWANT exactly by resetting TND with a call to RESET.
C  On the first step this matter is handled differently as explained above.
C
         IF (.NOT.FIRST) THEN
            CALL RESET(TWANT)
            CALL CHKFL(ASK,BADERR)
            IF (BADERR) GO TO 100
         END IF
      END IF
C
C  Process output, decide whether to take another step.
C
   40 CONTINUE
C
      IF (INTP) THEN
C
C  Interpolation is possible with this METHD.  The integration has
C  already reached T. If this is past TWANT, GOBACK is set .TRUE. and
C  the answers are obtained by interpolation.
C
         GOBACK = DIR*(T-TWANT) .GE. ZERO
         IF (GOBACK) THEN
            CALL INTRP(TWANT,'Both solution and derivative',NEQN,YGOT,
     &                 YPGOT,F,WORK,WORK(PRINTP),LNINTP)
            CALL CHKFL(ASK,BADERR)
            IF (BADERR) GO TO 100
            TGOT = TWANT
         END IF
      ELSE
C
C  Interpolation is not possible with this METHD, so output is obtained
C  by integrating to TWANT = TND.  Both YGOT(*) and YPGOT(*) are then 
C  already loaded with the solution at TWANT by CT.
C
         GOBACK = T .EQ. TWANT
         IF (GOBACK) TGOT = TWANT
      END IF
C
C  Updating of YMAX(*) is done here to account for the fact that when
C  interpolation is done, the integration goes past TGOT.  Note that YGOT(*)
C  is not defined until CT is called.  YMAX(*) was initialized at TSTRT
C  from values stored in WORK(*), so only needs to be updated for T
C  different from TSTRT.
      IF (T.NE.TSTRT) THEN
         DO 60 L = 1, NEQN
            YMAX(L) = MAX(YMAX(L),ABS(YGOT(L)))
   60    CONTINUE
      END IF
C
C  If done, go to the exit point.
      IF (GOBACK) GO TO 100
C
C  Take a step with CT in the direction of TND.  On exit, the solution is
C  advanced to TNOW.  The way CT is written, the approximate solution at
C  TNOW is available in both YGOT(*) and in WORK(*).  If output is obtained by
C  stepping to the end (TNOW = TWANT = TND), YGOT(*) can be returned directly.
C  If output is obtained by interpolation, the subroutine INTRP that does this
C  uses the values in WORK(*) for its computations and places the approximate
C  solution at TWANT in the array YGOT(*) for return to the calling program.
C  The approximate derivative is handled in the same way. TNOW is output from
C  CT and is actually a copy of T declared above in a common block.
C
      CALL CT(F,TNOW,YGOT,YPGOT,WORK,CFLAG)
      IER = CFLAG
C
C  A successful step by CT is indicated by CFLAG = 1 or = 2.
      IF (CFLAG.EQ.1) THEN
         GO TO 40
      ELSE IF (CFLAG.EQ.2) THEN
C
C  Supplement the warning message written in CT.
         NREC = 3
         WRITE (REC,'(A/A/A)')
     &' ** The last message was produced on a call to CT from UT.  ',
     &' ** In UT the appropriate action is to change to METHOD = 2,',
     &' ** or, if insufficient memory is available, to METHOD = 1. '
      ELSE IF (CFLAG.LE.6) THEN
         NREC = 1
         WRITE (REC,'(A)')
     &' ** The last message was produced on a call to CT from UT.'
      ELSE
         BADERR = .TRUE.
      END IF
      TGOT = T
C
C  Update YMAX(*) before the return.
      DO 80 L = 1, NEQN
         YMAX(L) = MAX(YMAX(L),ABS(YGOT(L)))
   80 CONTINUE
C
C  Exit point for UT.
C
  100 CONTINUE
C
      IF (BADERR) THEN
         IER = 911
         NREC = 4
         WRITE (REC,'(A/A/A/A)')
     &' ** An internal call by UT to a subroutine resulted in an  ',
     &' ** error that should not happen.  Check your program      ',
     &' ** carefully for array sizes, correct number of arguments,',
     &' ** type mismatches ... .'
      END IF
C
      TLAST = TGOT
C
C  All exits are done here after a call to RKMSG to report
C  what happened and set UFLAG.
C
      CALL RKMSG(IER,SRNAME,NREC,UFLAG)
C
      RETURN
      END
      SUBROUTINE STAT(TOTFCN,STPCST,WASTE,STPSOK,HNEXT)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code STAT and how it is used in
C  conjunction with the integrators CT and UT, you should study the 
C  document file rksuite.doc carefully before attempting to use the code. 
C  The following "Brief Reminder" is intended only to remind you of the 
C  meaning, type, and size requirements of the arguments.
C
C  STAT is called to obtain some details about the integration.
C
C  OUTPUT VARIABLES
C
C     TOTFCN    - INTEGER
C                 Total number of calls to F in the integration so far --
C                 a measure of the cost of the integration.
C     STPCST    - INTEGER
C                 Cost of a typical step with this METHOD measured in
C                 calls to F.
C     WASTE     - DOUBLE PRECISION
C                 The number of attempted steps that failed to meet the 
C                 local error requirement divided by the total number of 
C                 steps attempted so far in the integration.
C     STPSOK    - INTEGER
C                 The number of accepted steps.
C     HNEXT     - DOUBLE PRECISION
C                 The step size the integrator plans to use for the next step.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HNEXT, WASTE
      INTEGER           STPCST, STPSOK, TOTFCN
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='STAT')
      LOGICAL           ASK
      INTEGER           MINUS1, MINUS2
      PARAMETER         (ASK=.TRUE.,MINUS1=-1,MINUS2=-2)
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D+0)
C     .. Local Scalars ..
      INTEGER           FLAG, IER, NREC, STATE
C     .. External Subroutines ..
      EXTERNAL          RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE
C     .. Executable Statements ..
C
      IER = 1
      NREC = 0
C
C  Is it permissible to call STAT?
C
      CALL RKSIT(ASK,SRNAME,STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 20
      END IF
      IF (STATE.EQ.MINUS2) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A/A)')
     &' ** You have already made a call to STAT after a hard   ',
     &' ** failure was reported from the integrator. You cannot',
     &' ** call STAT again.'
         GO TO 20
      END IF
      CALL RKSIT(ASK,'CT',STATE)
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         IF (UTASK) THEN
            WRITE (REC,'(A)')
     &' ** You have not called UT, so you cannot use STAT.'
         ELSE
            WRITE (REC,'(A)')
     &' ** You have not called CT, so you cannot use STAT.'
         END IF
         GO TO 20
      END IF
C
C  Set flag so that the routine can only be called once after a hard 
C  failure from the integrator.
      IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2
C
      TOTFCN = SVNFCN + NFCN
      IF (ERASON) TOTFCN = TOTFCN + GNFCN
      STPCST = COST
      STPSOK = OKSTP
      IF (OKSTP.LE.1) THEN
         WASTE = ZERO
      ELSE
         WASTE = DBLE(FLSTP)/DBLE(FLSTP+OKSTP)
      END IF
      HNEXT = H
C
   20 CONTINUE
C
      CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
      RETURN
      END
      SUBROUTINE GLBERR(RMSERR,ERRMAX,TERRMX,WORK)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code GLBERR and how it is used in
C  conjunction with UT and CT to solve initial value problems, you should
C  study the document file rksuite.doc carefully before attempting to use 
C  the code.  The following "Brief Reminder" is intended only to remind you 
C  of the meaning, type, and size requirements of the arguments.
C
C  If ERRASS was set .TRUE. in the call to SETUP, then after any call to UT
C  or CT to advance the integration to TNOW or TWANT, the subroutine GLBERR
C  may be called to obtain an assessment of the true error of the integration.
C  At each step and for each solution component Y(L), a more accurate "true"
C  solution YT(L), an average magnitude "size(L)" of its size, and its error
C                abs(Y(L) - YT(L))/max("size(L)",THRES(L))
C  are computed.  The assessment returned in RMSERR(L) is the RMS (root-mean-
C  square) average of the error in the Lth solution component over all steps
C  of the integration from TSTART through TNOW.
C
C  OUTPUT VARIABLES
C
C     RMSERR(*) - DOUBLE PRECISION array of length NEQ
C                 RMSERR(L) approximates the RMS average of the true error 
C                 of the numerical solution for the Ith solution component,
C                 L = 1,2,...,NEQ.  The average is taken over all steps from
C                 TSTART to TNOW.
C     ERRMAX    - DOUBLE PRECISION
C                 The maximum (approximate) true error taken over all
C                 solution components and all steps from TSTART to TNOW.
C     TERRMX    - DOUBLE PRECISION
C                 First value of the independent variable where the
C                 (approximate) true error attains the maximum value ERRMAX.
C
C  WORKSPACE
C
C     WORK(*)   - DOUBLE PRECISION array as used in SETUP and UT or CT
C                 Do not alter the contents of this array.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ERRMAX, TERRMX
C     .. Array Arguments ..
      DOUBLE PRECISION  RMSERR(*), WORK(*)
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='GLBERR')
      LOGICAL           ASK
      PARAMETER         (ASK=.TRUE.)
      INTEGER           MINUS1, MINUS2
      PARAMETER         (MINUS1=-1,MINUS2=-2)
C     .. Local Scalars ..
      INTEGER           FLAG, IER, L, NREC, STATE
C     .. External Subroutines ..
      EXTERNAL          RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
C
      IER = 1
      NREC = 0
C
C  Is it permissible to call GLBERR?
C
      CALL RKSIT(ASK,SRNAME,STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 40
      END IF
      IF (STATE.EQ.MINUS2) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A/A)')
     &' ** You have already made a call to GLBERR after a hard ',
     &' ** failure was reported from the integrator. You cannot',
     &' ** call GLBERR again.'
         GO TO 40
      END IF
      CALL RKSIT(ASK,'CT',STATE)
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         IF (UTASK) THEN
            WRITE (REC,'(A)')
     &' ** You have not yet called UT, so you cannot call GLBERR.'
         ELSE
            WRITE (REC,'(A)')
     &' ** You have not yet called CT, so you cannot call GLBERR.'
         END IF
         GO TO 40
      END IF
C
C  Set flag so that the routine can only be called once after a hard 
C  failure from the integrator.
      IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2
C
C  Check that ERRASS was set properly for error assessment in SETUP.
C
      IF (.NOT.ERASON) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A/A)')
     &' ** No error assessment is available since you did not ',
     &' ** ask for it in your call to the routine SETUP.',
     &' ** Check your program carefully.'
         GO TO 40
      END IF
C
C  Check to see if the integrator has not actually taken a step.
C
      IF (OKSTP.EQ.0) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A/A)')
     &' ** The integrator has not actually taken any successful ',
     &' ** steps.  You cannot call GLBERR in this circumstance. ',
     &' ** Check your program carefully.'
         GO TO 40
      END IF
C
C  Compute RMS error and set output variables.
C
      ERRMAX = MAXERR
      TERRMX = LOCMAX
      DO 20 L = 1, NEQN
         RMSERR(L) = SQRT(WORK(PRZERR-1+L)/OKSTP)
   20 CONTINUE
C
   40 CONTINUE
C
      CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
      RETURN
      END
      SUBROUTINE CT(F,TNOW,YNOW,YPNOW,WORK,CFLAG)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code CT and how it is used in
C  conjunction with SETUP to solve initial value problems, you should study
C  the document file rksuite.doc carefully before attempting to use the code.
C  The following "Brief Reminder" is intended only to remind you of the
C  meaning, type, and size requirements of the arguments.
C
C  NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM:
C
C     F         - name of the subroutine for evaluating the differential
C                 equations.
C
C  The subroutine F must have the form
C
C  SUBROUTINE F(T,Y,YP)
C  DOUBLE PRECISION T,Y(*),YP(*)
C     Using the input values of the independent variable T and the solution
C     components Y(*), for each L = 1,2,...,NEQ evaluate the differential
C     equation for the derivative of the Lth solution component and place the
C     value in YP(L).  Do not alter the input values of T and Y(*).
C  RETURN
C  END
C
C  OUTPUT VARIABLES
C
C     TNOW      - DOUBLE PRECISION
C                 Current value of the independent variable.
C     YNOW(*)   - DOUBLE PRECISION array of length NEQ
C                 Approximation to the true solution at TNOW.
C     YPNOW(*)  - DOUBLE PRECISION array of length NEQ
C                 Approximation to the first derivative of the
C                 true solution at TNOW.
C
C  WORKSPACE
C
C     WORK(*)   - DOUBLE PRECISION array as used in SETUP
C                 Do not alter the contents of this array.
C
C  OUTPUT VARIABLE
C
C     CFLAG     - INTEGER
C
C                       SUCCESS.  A STEP WAS TAKEN TO TNOW.
C                 = 1 - Complete success.
C
C                       "SOFT" FAILURES
C                 = 2 - Warning:  You have obtained an answer by integrating
C                       to TEND (TNOW = TEND).  You have done this at least
C                       100 times, and monitoring of the computation reveals
C                       that this way of getting output has degraded the
C                       efficiency of the code. If you really need answers at
C                       so many specific points, it would be more efficient to
C                       get them with INTRP.  (If METHOD = 3, you would need
C                       to change METHOD and restart from TNOW, YNOW(*) by a
C                       call to SETUP.)  If you wish to continue as you are,
C                       you may.
C                 = 3 - Warning:  A considerable amount of work has been
C                       expended. To continue the integration, just call
C                       CT again.
C                 = 4 - Warning:  It appears that this problem is "stiff".
C                       You really should change to another code that is
C                       intended for such problems, but if you insist, you 
C                       can continue with CT by calling it again.
C
C                       "HARD" FAILURES
C                 = 5 - You are asking for too much accuracy. You cannot
C                       continue integrating this problem.
C                 = 6 - The global error assessment may not be reliable beyond
C                       the current point in the integration.  You cannot
C                       continue integrating this problem.
C
C                       "CATASTROPHIC" FAILURES
C                 = 911 - The nature of the catastrophe is reported on
C                         the standard output channel. Unless special
C                         provision was made in advance (see rksuite.doc),
C                         the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TNOW
      INTEGER           CFLAG
C     .. Array Arguments ..
      DOUBLE PRECISION  WORK(*), YNOW(*), YPNOW(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block for General Workspace Pointers ..
      INTEGER           PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      COMMON /RKCOM3/   PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      SAVE   /RKCOM3/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='CT')
      LOGICAL           ASK, TELL
      PARAMETER         (ASK=.TRUE.,TELL=.FALSE.)
      INTEGER           MINUS1, MINUS2
      PARAMETER         (MINUS1=-1,MINUS2=-2)
      INTEGER           MAXFCN
      PARAMETER         (MAXFCN=5000)
      DOUBLE PRECISION  ZERO, PT1, PT9, ONE, TWO, HUNDRD
      PARAMETER         (ZERO=0.0D+0,PT1=0.1D+0,PT9=0.9D+0,ONE=1.0D+0,
     &                  TWO=2.0D+0,HUNDRD=100.0D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  ALPHA, BETA, ERR, ERROLD, HAVG, HMIN, HTRY, TAU,
     &                  TEMP1, TEMP2, YPNORM
      INTEGER           IER, JFLSTP, L, NREC, NTEND, POINT, STATE, YNEW,
     &                  YPOLD
      LOGICAL           CHKEFF, FAILED, MAIN, PHASE1, PHASE2, PHASE3,
     &                  TOOMCH
C     .. External Subroutines ..
      EXTERNAL          RKMSG, RKSIT, STEP, STIFF, TRUERR
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MIN, SIGN
C     .. Save statement ..
      SAVE              JFLSTP, NTEND, ERROLD, HAVG, PHASE2, YNEW,
     &                  YPOLD, CHKEFF
C     .. Executable Statements ..
C
      IER = 1
      NREC = 0
C
C  Is it permissible to call CT?
C
      CALL RKSIT(ASK,'SETUP',STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 180
      END IF
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A)')
     &' ** You have not called SETUP, so you cannot use CT.'
         GO TO 180
      END IF
      IF (UTASK) THEN
         CALL RKSIT(ASK,'UT',STATE)
         IF (STATE.NE.MINUS2) THEN
            IER = 911
            NREC = 2
            WRITE (REC,'(A/A)')
     &' ** You have called CT after you specified in SETUP that ',
     &' ** you were going to use UT. This is not permitted.'
            UTASK = .FALSE.
            GO TO 180
         END IF
      END IF
      CALL RKSIT(ASK,SRNAME,STATE)
      IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A)')
     &' ** CT has already returned with a flag value of 5 or 6.',
     &' ** You cannot continue integrating this problem. You must ',
     &' ** call SETUP to start another problem.'
         GO TO 180
      END IF
C
      IF (FIRST) THEN
C
C  First call in an integration -- initialize everything.
C
         CHKEFF = .FALSE.
         NTEND = 0
         JFLSTP = 0
C
C  A scratch area of WORK(*) starting at PRSCR is used to hold two
C  arrays in this subroutine: the higher order approximate solution at
C  the end of a step and the approximate derivative of the solution at
C  the end of the last step. To make this clear, local pointers YNEW and
C  YPOLD are used.
         YNEW = PRSCR
         YPOLD = PRSCR
C
C  For this first step T was initialized to TSTRT in SETUP and the
C  starting values YSTART(*) were loaded into the area of WORK(*) reserved
C  for the current solution approximation starting at location PRY. The
C  derivative is now computed and stored in WORK(*) starting at PRYP.
C  Subsequently these arrays are copied to the output vectors YNOW(*)
C  and YPNOW(*).
         CALL F(T,WORK(PRY),WORK(PRYP))
         NFCN = NFCN + 1
         DO 20 L = 1, NEQN
            YNOW(L) = WORK(PRY-1+L)
            YPNOW(L) = WORK(PRYP-1+L)
   20    CONTINUE
C
C  Set dependent variables for error assessment.
         IF (ERASON) THEN
            DO 40 L = 1, NEQN
               WORK(PRZY-1+L) = YNOW(L)
               WORK(PRZYP-1+L) = YPNOW(L)
   40       CONTINUE
         END IF
C
C  The weights for the control of the error depend on the size of the
C  solution at the beginning and at the end of the step. On the first
C  step we do not have all this information. Whilst determining the
C  initial step size we initialize the weight vector to the larger of
C  abs(Y(L)) and the threshold for this component.
         DO 60 L = 1, NEQN
            WORK(PRWT-1+L) = MAX(ABS(YNOW(L)),WORK(PRTHRS-1+L))
   60    CONTINUE
C
C  If HSTRT is equal to zero, the code is to find an on-scale initial step
C  size H.  CT has an elaborate scheme of three phases for finding such an H,
C  and some preparations are made earlier.  In SETUP an upper bound is placed
C  on H that reflects the scale of the independent variable. When UTASK is
C  .TRUE., UT refines this bound using the first output point.  Here in CT
C  PHASE1 applies a rule of thumb based on the error control, the order of the
C  the formula, and the size of the initial slope to get a crude approximation
C  to an on-scale H.  PHASE2 may reduce H in the course of taking the first
C  step.  PHASE3 repeatedly adjusts H and retakes the first step until H is
C  on scale.
C
C  A guess for the magnitude of the first step size H can be provided to SETUP
C  as HSTART.  If it is too big or too small, it is ignored and the automatic
C  determination of an on-scale initial step size is activated.  If it is
C  acceptable, H is set to HSTART in SETUP.  Even when H is supplied to CT,
C  PHASE3 of the scheme for finding an on-scale initial step size is made
C  active so that the code can deal with a bad guess.
C
         PHASE1 = HSTRT .EQ. ZERO
         PHASE2 = PHASE1
         PHASE3 = .TRUE.
         IF (PHASE1) THEN
            H = ABS(H)
            YPNORM = ZERO
            DO 80 L = 1, NEQN
               IF (ABS(YNOW(L)).NE.ZERO) THEN
                  YPNORM = MAX(YPNORM,ABS(YPNOW(L))/WORK(PRWT-1+L))
               END IF
   80       CONTINUE
            TAU = TOLR**EXPON
            IF (H*YPNORM.GT.TAU) H = TAU/YPNORM
            HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND)))
            H = DIR*MAX(H,HMIN)
            PHASE1 = .FALSE.
         END IF
C
      ELSE
C
C Continuation call
C
         IF (LAST) THEN
            IER = 911
            NREC = 3
            WRITE (REC,'(A,D13.5,A/A/A)')
     &' ** You have already reached TEND ( = ', TND, ').',
     &' ** To integrate further with the same problem you must ',
     &' ** call the routine RESET with a new value of TEND.'
            GO TO 180
         END IF
      END IF
C
C  Begin computation of a step here.
C
      FAILED = .FALSE.
C
  100 CONTINUE
      H = SIGN(ABS(H),DIR)
C
C  Reduce the step size if necessary so that the code will not step
C  past TND.  "Look ahead" to prevent unnecessarily small step sizes.
      LAST = DIR*((T+H)-TND) .GE. ZERO
      IF (LAST) THEN
         H = TND - T
      ELSE IF (DIR*((T+TWO*H)-TND).GE.ZERO) THEN
         H = (TND-T)/TWO
      END IF
C
C  When the integrator is at T and attempts a step of H, the function
C  defining the differential equations will be evaluated at a number of
C  arguments between T and T+H.  If H is too small, these arguments cannot
C  be clearly distinguished in the precision available.
C
      HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(T+H)))
      IF (ABS(H).LT.HMIN) THEN
         IER = 5
         NREC = 3
         WRITE (REC,'(A/A,D13.5,A,D13.5/A)')
     &' ** In order to satisfy your error requirements CT would ',
     &' ** have to use a step size of ',H,' at TNOW = ',T,
     &' ** This is too small for the machine precision.'
         GO TO 180
      END IF
C
C  Monitor the impact of output on the efficiency of the integration.
C
      IF (CHKEFF) THEN
         NTEND = NTEND + 1
         IF (NTEND.GE.100 .AND. NTEND.GE.OKSTP/3) THEN
            IER = 2
            NREC = 6
            WRITE (REC,'(A/A/A/A/A/A)')
     &' ** More than 100 output points have been obtained by ',
     &' ** integrating to TEND.  They have been sufficiently close ',
     &' ** to one another that the efficiency of the integration has ',
     &' ** been degraded. It would probably be (much) more efficient ',
     &' ** to obtain output by interpolating with INTRP (after ',
     &' ** changing to METHOD=2 if you are using METHOD = 3).'
            NTEND = 0
            GO TO 180
         END IF
      END IF
C
C  Check for stiffness and for too much work.  Stiffness can be
C  checked only after a successful step.
C
      IF (.NOT.FAILED) THEN
C
C  Check for too much work.
         TOOMCH = NFCN .GT. MAXFCN
         IF (TOOMCH) THEN
            IER = 3
            NREC = 3
           WRITE (REC,'(A,I6,A/A/A)')
     &' ** Approximately ',MAXFCN,' function evaluations have been ',
     &' ** used to compute the solution since the integration ',
     &' ** started or since this message was last printed.'
C
C  After this warning message, NFCN is reset to permit the integration
C  to continue.  The total number of function evaluations in the primary
C  integration is SVNFCN + NFCN.
            SVNFCN = SVNFCN + NFCN
            NFCN = 0
         END IF
C
C  Check for stiffness.  NREC is passed on to STIFF because when
C  TOOMCH = .TRUE. and stiffness is diagnosed, the message about too
C  much work is augmented inside STIFF to explain that it is due to
C  stiffness.
         CALL STIFF(F,HAVG,JFLSTP,TOOMCH,MAXFCN,WORK,IER,NREC)
C
         IF (IER.NE.1) GO TO 180
      END IF
C
C  Take a step.  Whilst finding an on-scale H (PHASE2 = .TRUE.), the input
C  value of H might be reduced (repeatedly), but it will not be reduced
C  below HMIN.  The local error is estimated, a weight vector is formed,
C  and a weighted maximum norm, ERR, of the local error is returned.
C  The variable MAIN is input as .TRUE. to tell STEP that this is the
C  primary, or "main", integration.
C
C  H resides in the common block /RKCOM2/ which is used by both CT and STEP;
C  since it may be changed inside STEP, a local copy is made to ensure
C  portability of the code.
C
      MAIN = .TRUE.
      HTRY = H
      CALL STEP(F,NEQN,T,WORK(PRY),WORK(PRYP),WORK(PRSTGS),TOLR,HTRY,
     &          WORK(PRWT),WORK(YNEW),WORK(PRERST),ERR,MAIN,HMIN,
     &          WORK(PRTHRS),PHASE2)
      H = HTRY
C
C  Compare the norm of the local error to the tolerance.
C
      IF (ERR.GT.TOLR) THEN
C
C  Failed step.  Reduce the step size and try again.
C
C  First step:  Terminate PHASE3 of the search for an on-scale step size.
C               The step size is not on scale, so ERR may not be accurate;
C               reduce H by a fixed factor.  Failed attempts to take the
C               first step are not counted.
C  Later step:  Use ERR to compute an "optimal" reduction of H.  More than
C               one failure indicates a difficulty with the problem and an
C               ERR that may not be accurate, so reduce H by a fixed factor.
C
         IF (FIRST) THEN
            PHASE3 = .FALSE.
            ALPHA = RS1
         ELSE
            FLSTP = FLSTP + 1
            JFLSTP = JFLSTP + 1
            IF (FAILED) THEN
               ALPHA = RS1
            ELSE
               ALPHA = SAFETY*(TOLR/ERR)**EXPON
               ALPHA = MAX(ALPHA,RS1)
            END IF
         END IF
         H = ALPHA*H
         FAILED = .TRUE.
         GO TO 100
      END IF
C
C  Successful step.
C
C  Predict a step size appropriate for the next step.  After the first
C  step the prediction can be refined using an idea of H.A. Watts that
C  takes account of how well the prediction worked on the previous step.
      BETA = (ERR/TOLR)**EXPON
      IF (.NOT.FIRST) THEN
         TEMP1 = (ERR**EXPON)/H
         TEMP2 = (ERROLD**EXPON)/HOLD
         IF (TEMP1.LT.TEMP2*HUNDRD .AND. TEMP2.LT.TEMP1*HUNDRD) THEN
            BETA = BETA*(TEMP1/TEMP2)
         END IF
      END IF
      ALPHA = RS3
      IF (SAFETY.LT.BETA*ALPHA) ALPHA = SAFETY/BETA
C
C  On the first step a search is made for an on-scale step size.  PHASE2
C  of the scheme comes to an end here because a step size has been found
C  that is both successful and has a credible local error estimate. Except
C  in the special case that the first step is also the last, the step is
C  repeated in PHASE3 as long as an increase greater than RS2 appears
C  possible.  An increase as big as RS3 is permitted.  A step failure
C  terminates PHASE3.
C
      IF (FIRST) THEN
         PHASE2 = .FALSE.
         PHASE3 = PHASE3 .AND. .NOT. LAST .AND. (ALPHA.GT.RS2)
         IF (PHASE3) THEN
            H = ALPHA*H
            GO TO 100
         END IF
      END IF
C
C  After getting on scale, step size changes are more restricted.
      ALPHA = MIN(ALPHA,RS)
      IF (FAILED) ALPHA = MIN(ALPHA,ONE)
      ALPHA = MAX(ALPHA,RS1)
      HOLD = H
      H = ALPHA*H
C
C  For the diagnosis of stiffness, an average accepted step size, HAVG,
C  must be computed and SAVEd.
      IF (FIRST) THEN
         HAVG = HOLD
      ELSE
         HAVG = PT9*HAVG + PT1*HOLD
      END IF
C
      FIRST = .FALSE.
      ERROLD = ERR
      TOLD = T
C
C  Take care that T is set to precisely TND when the end of the
C  integration is reached.
      IF (LAST) THEN
         T = TND
      ELSE
         T = T + HOLD
      END IF
C
C  Increment counter on accepted steps.  Note that successful steps
C  that are repeated whilst getting on scale are not counted.
      OKSTP = OKSTP + 1
C
C  Advance the current solution and its derivative.  (Stored in WORK(*)
C  with the first location being PRY and PRYP, respectively.)  Update the
C  previous solution and its derivative.  (Stored in WORK(*) with the first
C  location being PRYOLD and YPOLD, respectively.)  Note that the previous
C  derivative will overwrite YNEW(*).
C
      DO 120 L = 1, NEQN
         WORK(PRYOLD-1+L) = WORK(PRY-1+L)
         WORK(PRY-1+L) = WORK(YNEW-1+L)
         WORK(YPOLD-1+L) = WORK(PRYP-1+L)
  120 CONTINUE
C
      IF (FSAL) THEN
C
C  When FSAL = .TRUE., YP(*) is the last stage of the step.
         POINT = PRSTGS + (LSTSTG-1)*NEQN
         DO 140 L = 1, NEQN
            WORK(PRYP-1+L) = WORK(POINT-1+L)
  140    CONTINUE
      ELSE
C
C  Call F to evaluate YP(*).
         CALL F(T,WORK(PRY),WORK(PRYP))
         NFCN = NFCN + 1
      END IF
C
C  If global error assessment is desired, advance the secondary
C  integration from TOLD to T.
C
      IF (ERASON) THEN
         CALL TRUERR(F,NEQN,WORK(PRY),TOLR,WORK(PRWT),WORK(PRZY),
     &               WORK(PRZYP),WORK(PRZERR),WORK(PRZYNU),WORK(PRZERS),
     &               WORK(PRZSTG),IER)
         IF (IER.EQ.6) THEN
C
C  The global error estimating procedure has broken down. Treat it as a
C  failed step. The solution and derivative are reset to their values at
C  the beginning of the step since the last valid error assessment refers
C  to them.
            OKSTP = OKSTP - 1
            ERASFL = .TRUE.
            LAST = .FALSE.
            T = TOLD
            H = HOLD
            DO 160 L = 1, NEQN
               WORK(PRY-1+L) = WORK(PRYOLD-1+L)
               WORK(PRYP-1+L) = WORK(YPOLD-1+L)
  160       CONTINUE
            IF (OKSTP.GT.1) THEN
               NREC = 2
               WRITE (REC,'(A/A,D13.5,A)')
     &' ** The global error assessment may not be reliable for T past ',
     &' ** TNOW = ',T,'.  The integration is being terminated.'
            ELSE
               NREC = 2
               WRITE (REC,'(A/A)')
     &' ** The global error assessment algorithm failed at the start',
     &' ** the integration.  The integration is being terminated.'
            END IF
            GO TO 180
         END IF
      END IF
C
C
C  Exit point for CT
C
  180 CONTINUE
C
C  Set the output variables and flag that interpolation is permitted
C
      IF (IER.LT.911) THEN
         TNOW = T
         LAST = TNOW .EQ. TND
         CHKEFF = LAST
         DO 200 L = 1, NEQN
            YNOW(L) = WORK(PRY-1+L)
            YPNOW(L) = WORK(PRYP-1+L)
  200    CONTINUE
         IF (IER.EQ.1) THEN
            STATE = MINUS2
            CALL RKSIT(TELL,'INTRP',STATE)
         END IF
      END IF
C
C  Call RKMSG to report what happened and set CFLAG.
C
      CALL RKMSG(IER,SRNAME,NREC,CFLAG)
C
      RETURN
      END
      SUBROUTINE INTRP(TWANT,REQEST,NWANT,YWANT,YPWANT,F,WORK,WRKINT,
     &                 LENINT)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code INTRP and how it is used in
C  conjunction with CT to solve initial value problems, you should study the
C  document file rksuite.doc carefully before attempting to use the code. The
C  following "Brief Reminder" is intended only to remind you of the meaning,
C  type, and size requirements of the arguments.
C
C  When integrating with METHOD = 1 or 2, answers may be obtained inexpensively
C  between steps by interpolation. INTRP is called after a step by CT from a
C  previous value of T, called TOLD below, to the current value of T to get
C  an answer at TWANT. You can specify any value of TWANT you wish, but
C  specifying a value outside the interval [TOLD,T] is likely to yield
C  answers with unsatisfactory accuracy.
C
C  INPUT VARIABLE
C
C     TWANT     - DOUBLE PRECISION
C                 The value of the independent variable where a solution
C                 is desired.
C
C  The interpolant is to be evaluated at TWANT to approximate the solution
C  and/or its first derivative there.  There are three cases:
C
C  INPUT VARIABLE
C
C     REQEST    - CHARACTER*(*)
C                 Only the first character of REQEST is significant.
C                 REQEST(1:1) = `S' or `s'- compute approximate `S'olution 
C                                           only.
C                             = `D' or `d'- compute approximate first 
C                                           `D'erivative of the solution only.
C                             = `B' or `b'- compute `B'oth approximate solution
C                                           and first derivative.
C                 Constraint: REQEST(1:1) must be `S',`s',`D',`d',`B' or `b'.
C
C  If you intend to interpolate at many points, you should arrange for the
C  the interesting components to be the first ones because the code
C  approximates only the first NWANT components.
C
C  INPUT VARIABLE
C
C     NWANT     - INTEGER
C                 Only the first NWANT components of the answer are to be
C                 computed.
C                 Constraint:  NEQ >= NWANT >= 1
C
C  OUTPUT VARIABLES
C
C     YWANT(*)  - DOUBLE PRECISION array of length NWANT
C                 Approximation to the first NWANT components of the true
C                 solution at TWANT when REQESTed.
C     YPWANT(*) - DOUBLE PRECISION array of length NWANT
C                 Approximation to the first NWANT components of the first
C                 derivative of the true solution at TWANT when REQESTed.
C
C  NAME DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM CALLING INTRP:
C
C     F         - name of the subroutine for evaluating the differential
C                 equations as provided to CT.
C
C  WORKSPACE
C
C     WORK(*)   - DOUBLE PRECISION array as used in SETUP and CT
C                 Do not alter the contents of this array.
C
C     WRKINT(*) - DOUBLE PRECISION array of length LENINT
C                 Do not alter the contents of this array.
C
C     LENINT    - INTEGER
C                 Length of WRKINT. If
C                 METHOD = 1, LENINT must be at least 1
C                        = 2, LENINT must be at least NEQ+MAX(NEQ,5*NWANT)
C                        = 3--CANNOT BE USED WITH THIS SUBROUTINE
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TWANT
      INTEGER           LENINT, NWANT
      CHARACTER*(*)     REQEST
C     .. Array Arguments ..
      DOUBLE PRECISION  WORK(*), WRKINT(*), YPWANT(*), YWANT(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block for General Workspace Pointers ..
      INTEGER           PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      COMMON /RKCOM3/   PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      SAVE   /RKCOM3/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='INTRP')
      LOGICAL           ASK
      INTEGER           PLUS1, MINUS1, MINUS2
      PARAMETER         (ASK=.TRUE.,PLUS1=1,MINUS1=-1,MINUS2=-2)
C     .. Local Scalars ..
      INTEGER           FLAG, ICHK, IER, NREC, NWNTSV, STARTP, STATE,
     &                  STATE1
      LOGICAL           ININTP, LEGALR
      CHARACTER         REQST1
C     .. External Subroutines ..
      EXTERNAL          EVALI, FORMI, RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         MAX
C     .. Save statement ..
      SAVE              NWNTSV, ININTP, STARTP
C     .. Data statements ..
      DATA              NWNTSV/MINUS1/
C     .. Executable Statements ..
C
      IER = 1
      NREC = 0
C
C  Is it permissible to call INTRP?
C
      CALL RKSIT(ASK,'CT',STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 20
      END IF
      IF (UTASK) THEN
         CALL RKSIT(ASK,'UT',STATE1)
         IF (STATE1.NE.MINUS2) THEN
            IER = 911
            NREC = 2
            WRITE (REC,'(A/A)')
     &' ** You have called INTRP after you specified to SETUP ',
     &' ** that you were going to use UT. This is not permitted.'
            GO TO 20
         END IF
      END IF
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A)')
     &' ** You have not called CT, so you cannot use INTRP.'
         GO TO 20
      END IF
      IF (STATE.GT.PLUS1) THEN
         IER = 911
         NREC = 2
         WRITE (REC,'(A/A)')
     &' ** CT has returned with a flag value greater than 1.',
     &' ** You cannot call INTRP in this circumstance.'
         GO TO 20
      END IF
C
C  Check input
C
      REQST1 = REQEST(1:1)
      LEGALR = REQST1 .EQ. 'S' .OR. REQST1 .EQ. 's' .OR.
     &         REQST1 .EQ. 'D' .OR. REQST1 .EQ. 'd' .OR.
     &         REQST1 .EQ. 'B' .OR. REQST1 .EQ. 'b'
      IF (.NOT.LEGALR) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A/A,A,A/A)')
     &' ** You have set the first character of ',
     &' ** REQEST to be ''',REQST1,'''. It must be one of ',
     &' ** ''S'',''s'',''D'',''d'',''B'' or ''b''.'
         GO TO 20
      END IF
C
      IF (NWANT.GT.NEQN) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A,I6,A/A,I6,A/A)')
     &' ** You have specified the value of NWANT to be ',NWANT,'. This',
     &' ** is greater than ',NEQN,', which is the number of equations ',
     &' ** in the system being integrated.'
         GO TO 20
      ELSE IF (NWANT.LT.1) THEN
         IER = 911
         NREC = 3
         WRITE (REC,'(A,I6,A/A/A)')
     &' ** You have specified the value of NWANT to be ',NWANT,', but ',
     &' ** this is less than 1. You cannot interpolate a zero or ',
     &' ** negative number of components.'
         GO TO 20
      END IF
C
      IF (METHD.EQ.1) THEN
         IF (LENINT.LT.1) THEN
            IER = 911
            NREC = 2
            WRITE (REC,'(A,I6,A/A)')
     &' ** You have specified LENINT to be ',LENINT,'.',
     &' ** This is too small. LENINT must be at least 1.'
            GO TO 20
         END IF
         STARTP = 1
      ELSE IF (METHD.EQ.2) THEN
         ICHK = NEQN + MAX(NEQN,5*NWANT)
         IF (LENINT.LT.ICHK) THEN
            IER = 911
            NREC = 3
            WRITE (REC,'(A,I6,A/A/A,I6,A)')
     &' ** You have specified LENINT to be ',LENINT,'. This is too',
     &' ** small. NINT must be at least NEQ + MAX(NEQ, 5*NWANT) ',
     &' ** which is ', ICHK,'.'
            GO TO 20
         END IF
         STARTP = NEQN + 1
      ELSE IF (METHD.EQ.3) THEN
         IER = 911
         NREC = 5
         WRITE (REC,'(A/A/A/A/A)')
     &' ** You have been using CT with METHOD = 3 to integrate your  ',
     &' ** equations. You have just called INTRP, but interpolation  ',
     &' ** is not available for this METHOD. Either use METHOD = 2,  ',
     &' ** for which interpolation is available, or use RESET to make',
     &' ** CT step exactly to the points where you want output.'
         GO TO 20
      END IF
C
C  Has the interpolant been initialised for this step?
C
      CALL RKSIT(ASK,SRNAME,STATE)
      ININTP = STATE .NE. MINUS2
C
C  Some initialization must be done before interpolation is possible.
C  To reduce the overhead, the interpolating polynomial is formed for
C  the first NWANT components.  In the unusual circumstance that NWANT
C  is changed while still interpolating within the span of the current
C  step, the scheme must be reinitialized to accomodate the additional
C  components.
C
      IF (.NOT.ININTP .OR. NWANT.NE.NWNTSV) THEN
C
C  At present the derivative of the solution at the previous step, YPOLD(*),
C  is stored in the scratch array area starting at PRSCR. In the case of
C  METHD = 1 we can overwrite the stages.
C
         IF (METHD.EQ.1) THEN
            CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD),
     &                 WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP,
     &                 WORK(PRSTGS),WORK(PRSTGS))
         ELSE
            CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD),
     &                 WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP,WRKINT,
     &                 WRKINT(STARTP))
         END IF
C
C  Set markers to show that interpolation has been initialized for
C  NWANT components.
         NWNTSV = NWANT
         ININTP = .TRUE.
      END IF
C
C  The actual evaluation of the interpolating polynomial and/or its first
C  derivative is done in EVALI.
C
      IF (METHD.EQ.1) THEN
         CALL EVALI(WORK(PRY),WORK(PRYP),WORK(PRSTGS),TWANT,REQEST,
     &              NWANT,YWANT,YPWANT)
      ELSE
         CALL EVALI(WORK(PRY),WORK(PRYP),WRKINT(STARTP),TWANT,REQEST,
     &              NWANT,YWANT,YPWANT)
      END IF
C
   20 CONTINUE
C
      CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
      RETURN
      END
      SUBROUTINE RESET(TENDNU)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  If you are not familiar with the code RESET and how it is used in
C  conjunction with CT to solve initial value problems, you should study the
C  document file rksuite.doc carefully before attempting to use the code. The
C  following "Brief Reminder" is intended only to remind you of the meaning,
C  type, and size requirements of the arguments.
C
C  The integration using CT proceeds from TSTART in the direction of TEND, and
C  is now at TNOW.  To reset TEND to a new value TENDNU, just call RESET with
C  TENDNU as the argument.  You must continue integrating in the same
C  direction, so the sign of (TENDNU - TNOW) must be the same as the sign of
C  (TEND - TSTART). To change direction you must restart by a call to SETUP.
C
C  INPUT VARIABLE
C
C     TENDNU    - DOUBLE PRECISION
C                 The new value of TEND.
C                 Constraint: TENDNU and TNOW must be clearly distinguishable
C                 in the precision used.  The sign of (TENDNU - TNOW) must be
C                 the same as the sign of (TEND - TSTART).
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TENDNU
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='RESET')
      LOGICAL           ASK
      INTEGER           MINUS1, MINUS2
      PARAMETER         (ASK=.TRUE.,MINUS1=-1,MINUS2=-2)
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  HMIN, TDIFF
      INTEGER           FLAG, IER, NREC, STATE, STATE1
C     .. External Subroutines ..
      EXTERNAL          RKMSG, RKSIT
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX
C     .. Executable Statements ..
      IER = 1
      NREC = 0
C
C  Is it permissible to call RESET?
C
      CALL RKSIT(ASK,'CT',STATE)
      IF (STATE.EQ.911) THEN
         IER = 912
         NREC = 1
         WRITE (REC,'(A)')
     &' ** A catastrophic error has already been detected elsewhere.'
         GO TO 20
      END IF
      IF (UTASK) THEN
         CALL RKSIT(ASK,'UT',STATE1)
         IF (STATE1.NE.MINUS2) THEN
            IER = 911
            NREC = 2
            WRITE (REC,'(A/A)')
     &' ** You have called RESET after you specified to SETUP that ',
     &' ** you were going to use UT. This is not permitted.'
            GO TO 20
         END IF
      END IF
      IF (STATE.EQ.MINUS1) THEN
         IER = 911
         NREC = 1
         WRITE (REC,'(A)')
     &' ** You have not called CT, so you cannot use RESET.'
         GO TO 20
      END IF
      IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
         IER = 911
         NREC = 2
         WRITE (REC,'(A,I1,A/A)')
     &' ** CT has returned with CFLAG =  ',STATE,'.',
     &' ** You cannot call RESET in this circumstance.'
         GO TO 20
      END IF
C
C  Check value of TENDNU
C
      IF (DIR.GT.ZERO .AND. TENDNU.LE.T) THEN
         IER = 911
         NREC = 4
         WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)')
     &' ** Integration is proceeding in the positive direction. The ',
     &' ** current value for the independent variable is ',T,
     &' ** and you have set TENDNU = ',TENDNU,'.  TENDNU must be ',
     &' ** greater than T.'
      ELSE IF (DIR.LT.ZERO .AND. TENDNU.GE.T) THEN
         IER = 911
         NREC = 4
         WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)')
     &' ** Integration is proceeding in the negative direction. The ',
     &' ** current value for the independent variable is ',T,
     &' ** and you have set TENDNU = ',TENDNU,'.  TENDNU must be ',
     &' ** less than T.'
      ELSE
         HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(TENDNU)))
         TDIFF = ABS(TENDNU-T)
         IF (TDIFF.LT.HMIN) THEN
            IER = 911
            NREC = 4
            WRITE (REC,'(A,D13.5,A/A,D13.5,A/A/A,D13.5,A)')
     &' ** The current value of the independent variable T is ',T,'.',
     &' ** The TENDNU you supplied has ABS(TENDNU-T) = ',TDIFF,'.',
     &' ** For the METHOD and the precision of the computer being ',
     &' ** used, this difference must be at least ',HMIN,'.'
         END IF
      END IF
      IF (IER.EQ.911) GO TO 20
C
      TND = TENDNU
      LAST = .FALSE.
C
   20 CONTINUE
C
      CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
      RETURN
      END
      SUBROUTINE MCONST(METHOD)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:   Sets machine-dependent global quantities
C
C  Common:    Initializes:    /RKCOM7/ OUTCH, MCHEPS, DWARF, RNDOFF,
C                                      SQRRMC, CUBRMC, TINY
C             Reads:          none
C             Alters:         none
C
C  Comments:
C  =========
C  OUTCH, MCHEPS, DWARF are pure environmental parameters with values
C  obtained from a call to ENVIRN. The other quantities set depend on
C  the environmental parameters, the implementation, and, possibly,
C  METHOD. At present the METHODs implemented in the RK suite do not
C  influence the values of these quantities.
C  OUTCH  - Standard output channel
C  MCHEPS - Largest positive number such that 1.0D0 + MCHEPS = 1.0D0
C  DWARF  - Smallest positive number
C  RNDOFF - 10 times MCHEPS
C  SQRRMC - Square root of MCHEPS
C  CUBRMC - Cube root of MCHEPS
C  TINY   - Square root of DWARF
C
C     .. Scalar Arguments ..
      INTEGER           METHOD
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Parameters ..
      DOUBLE PRECISION  TEN, THIRD
      PARAMETER         (TEN=10.0D+0,THIRD=1.0D+0/3.0D+0)
C     .. External Subroutines ..
      EXTERNAL          ENVIRN
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
C
      CALL ENVIRN(OUTCH,MCHEPS,DWARF)
C
      RNDOFF = TEN*MCHEPS
      SQRRMC = SQRT(MCHEPS)
      CUBRMC = MCHEPS**THIRD
      TINY = SQRT(DWARF)
C
      RETURN
      END
      SUBROUTINE ENVIRN(OUTCH,MCHEPS,DWARF)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  The RK suite requires some environmental parameters that are provided by
C  this subroutine.  The values provided with the distribution codes are those
C  appropriate to the IEEE standard.  They must be altered, if necessary, to
C  those appropriate to the computing system you are using before calling the 
C  codes of the suite.  
C
C        ================================================================
C        ================================================================
C        TO MAKE SURE THAT THESE MACHINE AND INSTALLATION DEPENDENT 
C        QUANTITIES ARE SPECIFIED PROPERLY, THE DISTRIBUTION VERSION 
C        WRITES A MESSAGE ABOUT THE MATTER TO THE STANDARD OUTPUT CHANNEL
C        AND TERMINATES THE RUN.  THE VALUES PROVIDED IN THE DISTRIBUTION
C        VERSION SHOULD BE ALTERED, IF NECESSARY, AND THE "WRITE" AND 
C        "STOP" STATEMENTS COMMENTED OUT.
C        ================================================================
C        ================================================================
C
C  OUTPUT VARIABLES 
C
C     OUTCH     - INTEGER
C                 Standard output channel
C     MCHEPS    - DOUBLE PRECISION
C                 MCHEPS is the largest positive number such that
C                 1.0D0 + MCHEPS = 1.0D0. 
C     DWARF     - DOUBLE PRECISION
C                 DWARF is the smallest positive number.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C     .. Scalar Arguments ..  
      INTEGER           OUTCH
      DOUBLE PRECISION  DWARF, MCHEPS
C     .. Executable Statements ..      
C
C  The following six statements are to be Commented out after verification that
C  the machine and installation dependent quantities are specified correctly.
C  If you pass copies of RKSUITE on to others, please give them the whole
C  distribution version of RKSUITE, and in particular, give them a version 
C  of ENVIRN that does not have the following six statements Commented out.
CC       WRITE(*,*) ' Before using RKSUITE, you must verify that the  '
CC       WRITE(*,*) ' machine- and installation-dependent quantities  '
CC       WRITE(*,*) ' specified in the subroutine ENVIRN are correct, '
CC       WRITE(*,*) ' and then Comment these WRITE statements and the '
CC       WRITE(*,*) ' STOP statement out of ENVIRN.                   '
CC       STOP
C
C  The following values are appropriate to IEEE arithmetic with the typical
C  standard output channel.
C
      OUTCH  = 6
      MCHEPS = 1.111D-16
      DWARF  = 2.23D-308 !we can actually do better but it is safe as is
C      
C------------------------------------------------------------------------------
C  If you have the routines D1MACH and I1MACH on your system, you could
C  replace the preceding statements by the following ones to obtain the 
C  appropriate machine dependent numbers. The routines D1MACH and I1MACH 
C  are public domain software.  They are available from NETLIB.
C      .. Scalar Arguments ..  
C      INTEGER           OUTCH
C      DOUBLE PRECISION  DWARF, MCHEPS
C      .. External Functions ..
C      INTEGER           I1MACH
C      DOUBLE PRECISION  D1MACH
C      .. Executable Statements ..
C
C      OUTCH = I1MACH(2)
C      MCHEPS = D1MACH(3)
C      DWARF = D1MACH(1)
C
C  If you have the NAG Fortran Library available on your system, you could 
C  replace the preceding statements by the following ones to obtain the 
C  appropriate machine dependent numbers.
C
C      .. Scalar Arguments ..  
C      INTEGER           OUTCH
C      DOUBLE PRECISION  DWARF, MCHEPS
C      .. External Functions ..
C      DOUBLE PRECISION  X02AJF, X02AMF
C      .. Executable Statements ..
C
C      CALL X04AAF(0,OUTCH)
C      MCHEPS = X02AJF()
C      DWARF = X02AMF()
C
C  If you have the IMSL MATH/LIBRARY available on your system, you could
C  replace the preceding statements by the following ones to obtain the
C  appropriate machine dependent numbers.
C
C      .. Scalar Arguments ..  
C      INTEGER           OUTCH
C      DOUBLE PRECISION  DWARF, MCHEPS
C      .. External Functions ..
C      DOUBLE PRECISION  DMACH
C      .. Executable Statements ..
C
C      CALL UMACH(2,OUTCH)
C      MCHEPS = DMACH(4)
C      DWARF = DMACH(1)
C------------------------------------------------------------------------------
C
      RETURN
      END
      SUBROUTINE CONST(METHOD,VECSTG,REQSTG,LINTPL)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
*************************************************
C
C  Purpose:   Set formula definitions and formula characteristics for
C             selected method. Return storage requirements for the
C             selected method.
C
C  Input:     METHOD
C  Output:    VECSTG, REQSTG, LINTPL
C
C  Common:    Initializes:    /RKCOM4/ A(*,*), B(*), C(*), BHAT(*), R(*),
C                                      E(*), PTR(*), NSTAGE, METHD, INTP, MINTP
C                             /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD,
C                                      TANANG, RS, RS1, RS2, RS3, RS4, ORDER,
C                                      LSTSTG, MAXTRY, NSEC, FSAL
C             Reads:          /RKCOM7/ RNDOFF
C             Alters:         none
C
C  Comments:
C  =========
C  Runge-Kutta formula pairs are described by a set of coefficients
C  and by the setting of a number of parameters that describe the
C  characteristics of the pair.  The input variable METHD indicates
C  which of the three pairs available is to be set. In the case of
C  METHD = 2 additional coefficients are defined that make interpolation
C  of the results possible and provide an additional error estimator.
C  VECSTG is the number of columns of workspace required to compute the
C  stages of a METHD. For interpolation purposes the routine returns via
C  COMMON the logical variable INTP indicating whether interpolation is
C  possible and via the call list:
C  REQSTG - whether the stages are required to form the
C           interpolant (set .FALSE. if INTP=.FALSE.)
C  LINTPL - the number of extra columns of storage required for use
C           with UT (set 0 if INTP=.FALSE.)
C
C  Quantities set in common blocks:
C  METHD - copy of METHOD
C  A, B, C, BHAT - coefficients of the selected method
C  R      - extra coefficents for interpolation with METHD = 2
C  E      - extra coefficients for additional local error estimate
C           with METHD = 2
C  PTR    - vector of pointers indicating how individual stages are to
C           be stored.  With it zero coefficients of the formulas can
C           be exploited to reduce the storage required
C  NSTAGE - number of stages for the specified METHD
C  INTP   - indicates whether there is an associated interpolant
C           (depending on the method, the user may have to supply
C           extra workspace)
C  MINTP  - the degree of the interpolating polynomial, if one exists
C  FSAL   - indicates whether the last stage of a step can be used as
C           the first stage of the following step
C  LSTSTG - pointer to location of last stage for use with FSAL=.TRUE.
C  ORDER  - the lower order of the pair of Runge-Kutta formulas that
C           constitute a METHD
C  TANANG, 
C  STBRAD - the stability region of the formula used to advance
C           the integration is approximated by a sector in the left half
C           complex plane.  TANANG is the tangent of the interior angle
C           of the sector and STBRAD is the radius of the sector.
C  COST   - cost of a successful step in function evaluations
C  MAXTRY - limit on the number of iterations in the stiffness check. As
C           set, no more than 24 function evaluations are made in the check.
C  NSEC   - each step of size H in the primary integration corresponds to
C           NSEC steps of size H/NSEC in the secondary integration when
C           global error assessment is done.
C  EXPON  - used to adjust the step size; this code implements an error
C           per step control for which EXPON = 1/(ORDER + 1).
C  SAFETY - quantity used in selecting the step size
C  TOOSML - quantity used to determine when a step size is too small for
C           the precision available
C  RS, RS1,
C  RS2, RS3,
C  RS4    - quantities used in determining the maximum and minimum change
C           change in step size (set independently of METHD)
C
C  Further comments on SAFETY:
C ============================
C  The code estimates the largest step size that will yield the specified
C  accuracy.  About half the time this step size would result in a local
C  error that is a little too large, and the step would be rejected.  For
C  this reason a SAFETY factor is used to reduce the "optimal" value to one
C  more likely to succeed.  Unfortunately, there is some art in choosing this
C  value. The more expensive a failed step, the smaller one is inclined to 
C  take this factor. However, a small factor means that if the prediction were
C  good, more accuracy than desired would be obtained and the behavior of the
C  error would then be irregular.  The more stringent the tolerance, the better
C  the prediction, except near limiting precision. Thus the general range of 
C  tolerances expected influences the choice of SAFETY.
C
C     .. Scalar Arguments ..
      INTEGER           LINTPL, METHOD, VECSTG
      LOGICAL           REQSTG
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Parameters ..
      DOUBLE PRECISION  ONE, ZERO, TWO, FIFTY, FIVEPC
      PARAMETER         (ONE=1.0D+0,ZERO=0.0D+0,TWO=2.0D+0,FIFTY=50.D+0,
     &                  FIVEPC=0.05D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  CDIFF, DIFF
      INTEGER           I, J
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, INT, MAX, MIN
C     .. Executable Statements ..
C
      METHD = METHOD
C
      GO TO (20,40,100) METHD
C
C  METHD = 1.
C    This pair is from "A 3(2) Pair of Runge-Kutta Formulas" by P. Bogacki
C    and L.F. Shampine, Appl. Math. Lett., 2, pp. 321-325, 1989.  The authors
C    are grateful to P. Bogacki for his assistance in implementing the pair.
C
   20 CONTINUE
      NSTAGE = 4
      FSAL = .TRUE.
      ORDER = 2
      TANANG = 8.9D0
      STBRAD = 2.3D0
      SAFETY = 0.8D0
      INTP = .TRUE.
      MINTP = 3
      REQSTG = .FALSE.
      LINTPL = 2
      NSEC = 3
C
      PTR(1) = 0
      PTR(2) = 1
      PTR(3) = 2
      PTR(4) = 3
C
      A(2,1) = 1.D0/2.D0
      A(3,1) = 0.D0
      A(3,2) = 3.D0/4.D0
      A(4,1) = 2.D0/9.D0
      A(4,2) = 1.D0/3.D0
      A(4,3) = 4.D0/9.D0
C
C  The coefficients BHAT(*) refer to the formula used to advance the
C  integration, here the one of order 3.  The coefficients B(*) refer
C  to the other formula, here the one of order 2. For this pair, BHAT(*)
C  is not needed since FSAL = .TRUE.
C
      B(1) = 7.D0/24.D0
      B(2) = 1.D0/4.D0
      B(3) = 1.D0/3.D0
      B(4) = 1.D0/8.D0
C
      C(1) = 0.D0
      C(2) = 1.D0/2.D0
      C(3) = 3.D0/4.D0
      C(4) = 1.D0
C
      GO TO 120
C
C  METHD = 2
C    This pair is from "An Efficient Runge-Kutta (4,5) Pair" by P. Bogacki
C    and L.F. Shampine, Rept. 89-20, Math. Dept., Southern Methodist
C    University, Dallas, Texas, USA, 1989.  The authors are grateful to
C    P. Bogacki for his assistance in implementing the pair.  Shampine and
C    Bogacki subsequently modified the formula to enhance the reliability of
C    the pair.  The original fourth order formula is used in an estimate of
C    the local error.  If the step fails, the computation is broken off.  If
C    the step is acceptable, the first evaluation of the next step is done,
C    i.e., the pair is implemented as FSAL and the local error of the step
C    is again estimated with a fourth order formula using the additional data.
C    The step must succeed with both estimators to be accepted.  When the
C    second estimate is formed, it is used for the subsequent adjustment of
C    the step size because it is of higher quality.  The two fourth order
C    formulas are well matched to leading order, and only exceptionally do
C    the estimators disagree -- problems with discontinuous coefficients are
C    handled more reliably by using two estimators as is global error
C    estimation.
C
   40 CONTINUE
      NSTAGE = 8
      FSAL = .TRUE.
      ORDER = 4
      TANANG = 5.2D0
      STBRAD = 3.9D0
      SAFETY = 0.8D0
      INTP = .TRUE.
      REQSTG = .TRUE.
      MINTP = 6
      LINTPL = 6
      NSEC = 2
C
      PTR(1) = 0
      PTR(2) = 1
      PTR(3) = 2
      PTR(4) = 3
      PTR(5) = 4
      PTR(6) = 5
      PTR(7) = 6
      PTR(8) = 7
C
      A(2,1) = 1.D0/6.D0
      A(3,1) = 2.D0/27.D0
      A(3,2) = 4.D0/27.D0
      A(4,1) = 183.D0/1372.D0
      A(4,2) = -162.D0/343.D0
      A(4,3) = 1053.D0/1372.D0
      A(5,1) = 68.D0/297.D0
      A(5,2) = -4.D0/11.D0
      A(5,3) = 42.D0/143.D0
      A(5,4) = 1960.D0/3861.D0
      A(6,1) = 597.D0/22528.D0
      A(6,2) = 81.D0/352.D0
      A(6,3) = 63099.D0/585728.D0
      A(6,4) = 58653.D0/366080.D0
      A(6,5) = 4617.D0/20480.D0
      A(7,1) = 174197.D0/959244.D0
      A(7,2) = -30942.D0/79937.D0
      A(7,3) = 8152137.D0/19744439.D0
      A(7,4) = 666106.D0/1039181.D0
      A(7,5) = -29421.D0/29068.D0
      A(7,6) = 482048.D0/414219.D0
      A(8,1) = 587.D0/8064.D0
      A(8,2) = 0.D0
      A(8,3) = 4440339.D0/15491840.D0
      A(8,4) = 24353.D0/124800.D0
      A(8,5) = 387.D0/44800.D0
      A(8,6) = 2152.D0/5985.D0
      A(8,7) = 7267.D0/94080.D0
C
C  The coefficients B(*) refer to the formula of order 4.
C
      B(1) = 2479.D0/34992.D0
      B(2) = 0.D0
      B(3) = 123.D0/416.D0
      B(4) = 612941.D0/3411720.D0
      B(5) = 43.D0/1440.D0
      B(6) = 2272.D0/6561.D0
      B(7) = 79937.D0/1113912.D0
      B(8) = 3293.D0/556956.D0
C
C  The coefficients E(*) refer to an estimate of the local error based on
C  the first formula of order 4.  It is the difference of the fifth order
C  result, here located in A(8,*), and the fourth order result.  By
C  construction both E(2) and E(7) are zero.
C
      E(1) = -3.D0/1280.D0
      E(2) = 0.D0
      E(3) = 6561.D0/632320.D0
      E(4) = -343.D0/20800.D0
      E(5) = 243.D0/12800.D0
      E(6) = -1.D0/95.D0
      E(7) = 0.D0
C
      C(1) = 0.D0
      C(2) = 1.D0/6.D0
      C(3) = 2.D0/9.D0
      C(4) = 3.D0/7.D0
      C(5) = 2.D0/3.D0
      C(6) = 3.D0/4.D0
      C(7) = 1.D0
      C(8) = 1.D0
C
C  To do interpolation with this pair, some extra stages have to be computed.
C  The following additional A(*,*) and C(*) coefficients are for this purpose.
C  In addition there is an array R(*,*) that plays a role for interpolation
C  analogous to that of BHAT(*) for the basic step.
C
      C(9) = 1.D0/2.D0
      C(10) = 5.D0/6.D0
      C(11) = 1.D0/9.D0
C
      A(9,1) = 455.D0/6144.D0
      A(10,1) = -837888343715.D0/13176988637184.D0
      A(11,1) = 98719073263.D0/1551965184000.D0
      A(9,2) = 0.D0
      A(10,2) = 30409415.D0/52955362.D0
      A(11,2) = 1307.D0/123552.D0
      A(9,3) = 10256301.D0/35409920.D0
      A(10,3) = -48321525963.D0/759168069632.D0
      A(11,3) = 4632066559387.D0/70181753241600.D0
      A(9,4) = 2307361.D0/17971200.D0
      A(10,4) = 8530738453321.D0/197654829557760.D0
      A(11,4) = 7828594302389.D0/382182512025600.D0
      A(9,5) = -387.D0/102400.D0
      A(10,5) = 1361640523001.D0/1626788720640.D0
      A(11,5) = 40763687.D0/11070259200.D0
      A(9,6) = 73.D0/5130.D0
      A(10,6) = -13143060689.D0/38604458898.D0
      A(11,6) = 34872732407.D0/224610586200.D0
      A(9,7) = -7267.D0/215040.D0
      A(10,7) = 18700221969.D0/379584034816.D0
      A(11,7) = -2561897.D0/30105600.D0
      A(9,8) = 1.D0/32.D0
      A(10,8) = -5831595.D0/847285792.D0
      A(11,8) = 1.D0/10.D0
      A(10,9) = -5183640.D0/26477681.D0
      A(11,9) = -1.D0/10.D0
      A(11,10) = -1403317093.D0/11371610250.D0
C
      DO 60 I = 1, 11
         R(I,1) = 0.D0
   60 CONTINUE
      DO 80 I = 1, 6
         R(2,I) = 0.D0
   80 CONTINUE
      R(1,6) = -12134338393.D0/1050809760.D0
      R(1,5) = -1620741229.D0/50038560.D0
      R(1,4) = -2048058893.D0/59875200.D0
      R(1,3) = -87098480009.D0/5254048800.D0
      R(1,2) = -11513270273.D0/3502699200.D0
C
      R(3,6) = -33197340367.D0/1218433216.D0
      R(3,5) = -539868024987.D0/6092166080.D0
      R(3,4) = -39991188681.D0/374902528.D0
      R(3,3) = -69509738227.D0/1218433216.D0
      R(3,2) = -29327744613.D0/2436866432.D0
C
      R(4,6) = -284800997201.D0/19905339168.D0
      R(4,5) = -7896875450471.D0/165877826400.D0
      R(4,4) = -333945812879.D0/5671036800.D0
      R(4,3) = -16209923456237.D0/497633479200.D0
      R(4,2) = -2382590741699.D0/331755652800.D0
C
      R(5,6) = -540919.D0/741312.D0
      R(5,5) = -103626067.D0/43243200.D0
      R(5,4) = -633779.D0/211200.D0
      R(5,3) = -32406787.D0/18532800.D0
      R(5,2) = -36591193.D0/86486400.D0
C
      R(6,6) = 7157998304.D0/374350977.D0
      R(6,5) = 30405842464.D0/623918295.D0
      R(6,4) = 183022264.D0/5332635.D0
      R(6,3) = -3357024032.D0/1871754885.D0
      R(6,2) = -611586736.D0/89131185.D0
C
      R(7,6) = -138073.D0/9408.D0
      R(7,5) = -719433.D0/15680.D0
      R(7,4) = -1620541.D0/31360.D0
      R(7,3) = -385151.D0/15680.D0
      R(7,2) = -65403.D0/15680.D0
C
      R(8,6) = 1245.D0/64.D0
      R(8,5) = 3991.D0/64.D0
      R(8,4) = 4715.D0/64.D0
      R(8,3) = 2501.D0/64.D0
      R(8,2) = 149.D0/16.D0
      R(8,1) = 1.D0
C
      R(9,6) = 55.D0/3.D0
      R(9,5) = 71.D0
      R(9,4) = 103.D0
      R(9,3) = 199.D0/3.D0
      R(9,2) = 16.0D0
C
      R(10,6) = -1774004627.D0/75810735.D0
      R(10,5) = -1774004627.D0/25270245.D0
      R(10,4) = -26477681.D0/359975.D0
      R(10,3) = -11411880511.D0/379053675.D0
      R(10,2) = -423642896.D0/126351225.D0
C
      R(11,6) = 35.D0
      R(11,5) = 105.D0
      R(11,4) = 117.D0
      R(11,3) = 59.D0
      R(11,2) = 12.D0
C
      GO TO 120
C
C  METHD = 3
C    This pair is from "High Order Embedded Runge-Kutta Formulae" by P.J.
C    Prince and J.R. Dormand, J. Comp. Appl. Math.,7, pp. 67-75, 1981.  The
C    authors are grateful to P. Prince and J. Dormand for their assistance in
C    implementing the pair.
C
  100 CONTINUE
      NSTAGE = 13
      FSAL = .FALSE.
      ORDER = 7
      TANANG = 11.0D0
      STBRAD = 5.2D0
      SAFETY = 0.8D0
      INTP = .FALSE.
      REQSTG = .FALSE.
      MINTP = 0
      LINTPL = 0
      NSEC = 2
C
      PTR(1) = 0
      PTR(2) = 1
      PTR(3) = 2
      PTR(4) = 1
      PTR(5) = 3
      PTR(6) = 2
      PTR(7) = 4
      PTR(8) = 5
      PTR(9) = 6
      PTR(10) = 7
      PTR(11) = 8
      PTR(12) = 9
      PTR(13) = 1
C
      A(2,1) = 5.55555555555555555555555555556D-2
      A(3,1) = 2.08333333333333333333333333333D-2
      A(3,2) = 6.25D-2
      A(4,1) = 3.125D-2
      A(4,2) = 0.D0
      A(4,3) = 9.375D-2
      A(5,1) = 3.125D-1
      A(5,2) = 0.D0
      A(5,3) = -1.171875D0
      A(5,4) = 1.171875D0
      A(6,1) = 3.75D-2
      A(6,2) = 0.D0
      A(6,3) = 0.D0
      A(6,4) = 1.875D-1
      A(6,5) = 1.5D-1
      A(7,1) = 4.79101371111111111111111111111D-2
      A(7,2) = 0.D0
      A(7,3) = 0.0D0
      A(7,4) = 1.12248712777777777777777777778D-1
      A(7,5) = -2.55056737777777777777777777778D-2
      A(7,6) = 1.28468238888888888888888888889D-2
      A(8,1) = 1.6917989787292281181431107136D-2
      A(8,2) = 0.D0
      A(8,3) = 0.D0
      A(8,4) = 3.87848278486043169526545744159D-1
      A(8,5) = 3.59773698515003278967008896348D-2
      A(8,6) = 1.96970214215666060156715256072D-1
      A(8,7) = -1.72713852340501838761392997002D-1
      A(9,1) = 6.90957533591923006485645489846D-2
      A(9,2) = 0.D0
      A(9,3) = 0.D0
      A(9,4) = -6.34247976728854151882807874972D-1
      A(9,5) = -1.61197575224604080366876923982D-1
      A(9,6) = 1.38650309458825255419866950133D-1
      A(9,7) = 9.4092861403575626972423968413D-1
      A(9,8) = 2.11636326481943981855372117132D-1
      A(10,1) = 1.83556996839045385489806023537D-1
      A(10,2) = 0.D0
      A(10,3) = 0.D0
      A(10,4) = -2.46876808431559245274431575997D0
      A(10,5) = -2.91286887816300456388002572804D-1
      A(10,6) = -2.6473020233117375688439799466D-2
      A(10,7) = 2.84783876419280044916451825422D0
      A(10,8) = 2.81387331469849792539403641827D-1
      A(10,9) = 1.23744899863314657627030212664D-1
      A(11,1) = -1.21542481739588805916051052503D0
      A(11,2) = 0.D0
      A(11,3) = 0.D0
      A(11,4) = 1.66726086659457724322804132886D1
      A(11,5) = 9.15741828416817960595718650451D-1
      A(11,6) = -6.05660580435747094755450554309D0
      A(11,7) = -1.60035735941561781118417064101D1
      A(11,8) = 1.4849303086297662557545391898D1
      A(11,9) = -1.33715757352898493182930413962D1
      A(11,10) = 5.13418264817963793317325361166D0
      A(12,1) = 2.58860916438264283815730932232D-1
      A(12,2) = 0.D0
      A(12,3) = 0.D0
      A(12,4) = -4.77448578548920511231011750971D0
      A(12,5) = -4.3509301377703250944070041181D-1
      A(12,6) = -3.04948333207224150956051286631D0
      A(12,7) = 5.57792003993609911742367663447D0
      A(12,8) = 6.15583158986104009733868912669D0
      A(12,9) = -5.06210458673693837007740643391D0
      A(12,10) = 2.19392617318067906127491429047D0
      A(12,11) = 1.34627998659334941535726237887D-1
      A(13,1) = 8.22427599626507477963168204773D-1
      A(13,2) = 0.D0
      A(13,3) = 0.D0
      A(13,4) = -1.16586732572776642839765530355D1
      A(13,5) = -7.57622116690936195881116154088D-1
      A(13,6) = 7.13973588159581527978269282765D-1
      A(13,7) = 1.20757749868900567395661704486D1
      A(13,8) = -2.12765911392040265639082085897D0
      A(13,9) = 1.99016620704895541832807169835D0
      A(13,10) = -2.34286471544040292660294691857D-1
      A(13,11) = 1.7589857770794226507310510589D-1
      A(13,12) = 0.D0
C
C  The coefficients BHAT(*) refer to the formula used to advance the
C  integration, here the one of order 8.  The coefficients B(*) refer
C  to the other formula, here the one of order 7.
C
      BHAT(1) = 4.17474911415302462220859284685D-2
      BHAT(2) = 0.D0
      BHAT(3) = 0.D0
      BHAT(4) = 0.D0
      BHAT(5) = 0.D0
      BHAT(6) = -5.54523286112393089615218946547D-2
      BHAT(7) = 2.39312807201180097046747354249D-1
      BHAT(8) = 7.0351066940344302305804641089D-1
      BHAT(9) = -7.59759613814460929884487677085D-1
      BHAT(10) = 6.60563030922286341461378594838D-1
      BHAT(11) = 1.58187482510123335529614838601D-1
      BHAT(12) = -2.38109538752862804471863555306D-1
      BHAT(13) = 2.5D-1
C
      B(1) = 2.9553213676353496981964883112D-2
      B(2) = 0.D0
      B(3) = 0.D0
      B(4) = 0.D0
      B(5) = 0.D0
      B(6) = -8.28606276487797039766805612689D-1
      B(7) = 3.11240900051118327929913751627D-1
      B(8) = 2.46734519059988698196468570407D0
      B(9) = -2.54694165184190873912738007542D0
      B(10) = 1.44354858367677524030187495069D0
      B(11) = 7.94155958811272872713019541622D-2
      B(12) = 4.44444444444444444444444444445D-2
      B(13) = 0.D0
C
      C(1) = 0.D0
      C(2) = 5.55555555555555555555555555556D-2
      C(3) = 8.33333333333333333333333333334D-2
      C(4) = 1.25D-1
      C(5) = 3.125D-1
      C(6) = 3.75D-1
      C(7) = 1.475D-1
      C(8) = 4.65D-1
      C(9) = 5.64865451382259575398358501426D-1
      C(10) = 6.5D-1
      C(11) = 9.24656277640504446745013574318D-1
      C(12) = 1.D0
      C(13) = C(12)
C
      GO TO 120
C
C  The definitions of all pairs come here for the calculation of
C  LSTSTG, RS1, RS2, RS3, RS4, COST, MAXTRY, EXPON, TOOSML, and VECSTG.
C
  120 CONTINUE
      LSTSTG = PTR(NSTAGE)
      IF (FSAL) THEN
         COST = DBLE(NSTAGE-1)
      ELSE
         COST = DBLE(NSTAGE)
      END IF
C
C  MAXTRY - limit on the number of iterations of a computation made in
C  diagnosing stiffness.  There are at most Q = 3 function calls per
C  iteration. MAXTRY is determined so that  Q*MAXTRY <= 5% of the cost of
C  50 steps and 1 <= MAXTRY <= 8. This limits the number of calls to FCN
C  in each diagnosis of stiffness to 24 calls.
C
      MAXTRY = MIN(8,MAX(1,INT(FIVEPC*COST*FIFTY)))
C
      EXPON = ONE/(ORDER+ONE)
C
C     In calculating CDIFF it is assumed that there will be a non-zero
C     difference |C(I) - C(J)| less than one. If C(I) = C(J) for any I not
C     equal to J, they should be made precisely equal by assignment.
C
      CDIFF = ONE
      DO 160 I = 1, NSTAGE - 1
         DO 140 J = I + 1, NSTAGE
            DIFF = ABS(C(I)-C(J))
            IF (DIFF.NE.ZERO) CDIFF = MIN(CDIFF,DIFF)
  140    CONTINUE
  160 CONTINUE
      TOOSML = RNDOFF/CDIFF
C
C  Determine the number of columns needed in STAGES(1:NEQ,*) (this will be
C  at most NSTAGE-1 since the first stage is held in a separate array).
C  The PTR array contains the column positions of the stages.
C
      VECSTG = 0
      DO 180 I = 2, NSTAGE
         VECSTG = MAX(PTR(I),VECSTG)
  180 CONTINUE
C
      RS = TWO
      RS1 = ONE/RS
      RS2 = RS**2
      RS3 = RS**3
      RS4 = ONE/RS3
C
      RETURN
      END
      SUBROUTINE FORMI(F,NEQ,NWANT,Y,YP,YOLD,YPOLD,STAGES,CALSTG,
     &                 XSTAGE,P)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:    Forms an interpolating polynomial for use with
C              METHDs 1 or 2.
C
C  Input:      NEQ, NWANT, T, Y(*), YP(*), HOLD, YOLD(*), YPOLD(*),
C              STAGES(NEQ,*), CALSTG
C  Output:     P(*), XSTAGE(NEQ)
C  External:   F
C
C  Common:     Initializes:    none
C              Reads:          /RKCOM4/ A(*,*), C(*), R(*), METHD, MINTP
C                              /RKCOM2/ T, TOLD, HOLD
C              Alters:         /RKCOM2/ NFCN
C
C  Comments:
C  =========
C  The integration has reached T with a step HOLD from TOLD = T-HOLD.
C  Y(*),YP(*) and YOLD(*),YPOLD(*) approximate the solution and its
C  derivative at T and TOLD respectively.  STAGES(NEQ,*) holds the stages
C  computed in taking this step. In the case of METHD = 2 it is necessary 
C  to compute some more stages in this subroutine. CALSTG indicates whether
C  or not the extra stages need to be computed. A(*,*) and C(*) are used in 
C  computing these stages. The extra stages are stored in STAGES(NEQ,*) and 
C  XSTAGE(*).  The coefficients of the interpolating polynomials for the first
C  NWANT components of the solution are returned in the array P(*). The 
C  polynomial is of degree MINTP = 3 for METHD = 1 and of degree MINTP = 6 
C  for METHD = 2. The vector R(*) is used for workspace when METHD = 2.
C
C     .. Scalar Arguments ..
      INTEGER           NEQ, NWANT
      LOGICAL           CALSTG
C     .. Array Arguments ..
      DOUBLE PRECISION  P(*), STAGES(NEQ,*), XSTAGE(*), Y(*), YOLD(*),
     &                  YP(*), YPOLD(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Local Scalars ..
      DOUBLE PRECISION  D1, D2, D3, D4, HYP, HYPOLD
      INTEGER           I, J, K, L
C     .. Executable Statements ..
C
      IF (METHD.EQ.1) THEN
C
C  METHD = 1.  Use the cubic Hermite interpolant that is is fully
C  specified by the values and slopes at the two ends of the step.
C
         DO 20 L = 1, NWANT
            D1 = Y(L) - YOLD(L)
            HYP = HOLD*YP(L)
            HYPOLD = HOLD*YPOLD(L)
            D2 = HYP - D1
            D3 = D1 - HYPOLD
            D4 = D2 - D3
            P(L) = D2 + D4
            P(NWANT+L) = D4
   20    CONTINUE
C
      ELSE
C
C  METHD = 2.
C
         IF (CALSTG) THEN
C
C  Compute the extra stages needed for interpolation using the facts that
C       1. Stage 1 is YPOLD(*).
C       2. Stage i (i>1) is stored in STAGES(1:NEQ,i).
C       3. This pair is FSAL, i.e. STAGES(1:NEQ,7)=YP(1:NEQ), which frees
C          up STAGES(1:NEQ,7) for use by stage 9.
C       4. XSTAGE(1:NEQ) is used for stage 10.
C       5. The coefficient of stage 2 in the interpolant is always 0, so
C          STAGES(1:NEQ,1) is used for stage 11.
C  The vector P(1:NEQ) is used as workspace for computing the stages.
C
            DO 180 I = 9, 11
               DO 40 L = 1, NEQ
                  P(L) = A(I,1)*YPOLD(L)
   40          CONTINUE
               DO 140 J = 2, I - 1
                  IF (J.LE.7) THEN
                     DO 60 L = 1, NEQ
                        P(L) = P(L) + A(I,J)*STAGES(L,J-1)
   60                CONTINUE
                  ELSE IF (J.EQ.8) THEN
                     DO 80 L = 1, NEQ
                        P(L) = P(L) + A(I,J)*YP(L)
   80                CONTINUE
                  ELSE IF (J.EQ.9) THEN
                     DO 100 L = 1, NEQ
                        P(L) = P(L) + A(I,J)*STAGES(L,7)
  100                CONTINUE
                  ELSE IF (J.EQ.10) THEN
                     DO 120 L = 1, NEQ
                        P(L) = P(L) + A(I,J)*XSTAGE(L)
  120                CONTINUE
                  END IF
  140          CONTINUE
               DO 160 L = 1, NEQ
                  P(L) = YOLD(L) + HOLD*P(L)
  160          CONTINUE
               IF (I.EQ.9) THEN
                  CALL F(TOLD+C(I)*HOLD,P,STAGES(1,7))
                  NFCN = NFCN + 1
               ELSE IF (I.EQ.10) THEN
                  CALL F(TOLD+C(I)*HOLD,P,XSTAGE)
                  NFCN = NFCN + 1
               ELSE
                  CALL F(TOLD+C(I)*HOLD,P,STAGES(1,1))
                  NFCN = NFCN + 1
               END IF
  180       CONTINUE
         END IF
C
C  Form the coefficients of the interpolating polynomial in its shifted
C  and scaled form.  The transformation from the form in which the
C  polynomial is derived can be somewhat ill-conditioned.  The terms 
C  are grouped so as to minimize the errors of the transformation.
C
C  Coefficient of SIGMA**6
         K = 4*NWANT
         DO 200 L = 1, NWANT
            P(K+L) = R(5,6)*STAGES(L,4) +
     &               ((R(10,6)*XSTAGE(L)+R(8,6)*YP(L))+
     &               (R(7,6)*STAGES(L,6)+R(6,6)*STAGES(L,5))) +
     &               ((R(4,6)*STAGES(L,3)+R(9,6)*STAGES(L,7))+
     &               (R(3,6)*STAGES(L,2)+R(11,6)*STAGES(L,1))+
     &               R(1,6)*YPOLD(L))
  200    CONTINUE
C
C  Coefficient of SIGMA**5
         K = 3*NWANT
         DO 220 L = 1, NWANT
            P(K+L) = (R(10,5)*XSTAGE(L)+R(9,5)*STAGES(L,7)) +
     &               ((R(7,5)*STAGES(L,6)+R(6,5)*STAGES(L,5))+
     &               R(5,5)*STAGES(L,4)) + ((R(4,5)*STAGES(L,3)+
     &               R(8,5)*YP(L))+(R(3,5)*STAGES(L,2)+R(11,5)*
     &               STAGES(L,1))+R(1,5)*YPOLD(L))
  220    CONTINUE
C
C  Coefficient of SIGMA**4
         K = 2*NWANT
         DO 240 L = 1, NWANT
            P(K+L) = ((R(4,4)*STAGES(L,3)+R(8,4)*YP(L))+
     &               (R(7,4)*STAGES(L,6)+R(6,4)*STAGES(L,5))+
     &               R(5,4)*STAGES(L,4)) + ((R(10,4)*XSTAGE(L)+
     &               R(9,4)*STAGES(L,7))+(R(3,4)*STAGES(L,2)+
     &               R(11,4)*STAGES(L,1))+R(1,4)*YPOLD(L))
  240    CONTINUE
C
C  Coefficient of SIGMA**3
         K = NWANT
         DO 260 L = 1, NWANT
            P(K+L) = R(5,3)*STAGES(L,4) + R(6,3)*STAGES(L,5) +
     &               ((R(3,3)*STAGES(L,2)+R(9,3)*STAGES(L,7))+
     &               (R(10,3)*XSTAGE(L)+R(8,3)*YP(L))+R(1,3)*
     &               YPOLD(L))+((R(4,3)*STAGES(L,3)+R(11,3)*
     &               STAGES(L,1))+R(7,3)*STAGES(L,6))
  260    CONTINUE
C
C  Coefficient of SIGMA**2
C
         DO 280 L = 1, NWANT
            P(L) = R(5,2)*STAGES(L,4) + ((R(6,2)*STAGES(L,5)+
     &             R(8,2)*YP(L))+R(1,2)*YPOLD(L)) +
     &             ((R(3,2)*STAGES(L,2)+R(9,2)*STAGES(L,7))+
     &             R(10,2)*XSTAGE(L)) + ((R(4,2)*STAGES(L,3)+
     &             R(11,2)*STAGES(L,1))+R(7,2)*STAGES(L,6))
  280    CONTINUE
C
C  Scale all the coefficients by the step size.
C
         DO 300 L = 1, NWANT*(MINTP-1)
            P(L) = HOLD*P(L)
  300    CONTINUE
C
      END IF
C
      RETURN
      END
      SUBROUTINE EVALI(Y,YP,P,TWANT,REQEST,NWANT,YWANT,YPWANT)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
*************************************************
C
C  Purpose:    Evaluation of an interpolating polynomial and/or its
C              first derivative.
C
C  Input:      Y(*), YP(*), P(NWANT,*), TWANT, REQEST, NWANT
C  Output:     YWANT(*), YPWANT(*)
C
C  Common:     Initializes:    none
C              Reads:          /RKCOM2/ HOLD, T
C                              /RKCOM4/ MINTP
C              Alters:         none
C
C  Comments:
C  =========
C  The interpolant is evaluated at TWANT to approximate the solution,
C  YWANT, and/or its first derivative there, YPWANT. Only the first
C  NWANT components of the answer are computed. There are three cases
C  that are indicated by the first character of REQEST:
C    REQEST(1:1) = `S' or `s'- compute approximate `S'olution only.
C                = `D' or `d'- compute approximate first `D'erivative
C                              of the solution only.
C                = `B' or `b'- compute `B'oth approximate solution and
C                              first derivative.
C  The coefficents of the polynomial are contained in Y(*), YP(*) and
C  P(NWANT,*).
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TWANT
      INTEGER           NWANT
      CHARACTER*(*)     REQEST
C     .. Array Arguments ..
      DOUBLE PRECISION  P(NWANT,*), Y(*), YP(*), YPWANT(*), YWANT(*)
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Local Scalars ..
      DOUBLE PRECISION  SIGMA
      INTEGER           K, L
      CHARACTER         REQST1
C     .. Executable Statements ..
C
C  Evaluate the interpolating polynomial of degree MINTP in terms of the
C  shifted and scaled independent variable SIGMA.
C
      SIGMA = (TWANT-T)/HOLD
C
      REQST1 = REQEST(1:1)
      IF (REQST1.EQ.'S' .OR. REQST1.EQ.'s' .OR. 
     &    REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN
C
         DO 20 L = 1, NWANT
            YWANT(L) = P(L,MINTP-1)*SIGMA
   20    CONTINUE
         DO 60 K = MINTP - 2, 1, -1
            DO 40 L = 1, NWANT
               YWANT(L) = (YWANT(L)+P(L,K))*SIGMA
   40       CONTINUE
   60    CONTINUE
         DO 80 L = 1, NWANT
            YWANT(L) = (YWANT(L)+HOLD*YP(L))*SIGMA + Y(L)
   80    CONTINUE
      END IF
C
C  Evaluate the derivative of the interpolating polynomial.
C
      IF (REQST1.EQ.'D' .OR. REQST1.EQ.'d' .OR. 
     &    REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN
C
C  The derivative of the interpolating polynomial with respect to TWANT 
C  is the derivative with respect to S divided by HOLD.
C
         DO 100 L = 1, NWANT
            YPWANT(L) = MINTP*P(L,MINTP-1)*SIGMA
  100    CONTINUE
         DO 140 K = MINTP - 1, 2, -1
            DO 120 L = 1, NWANT
               YPWANT(L) = (YPWANT(L)+K*P(L,K-1))*SIGMA
  120       CONTINUE
  140    CONTINUE
         DO 160 L = 1, NWANT
            YPWANT(L) = (YPWANT(L)+HOLD*YP(L))/HOLD
  160    CONTINUE
      END IF
C
      RETURN
      END
      SUBROUTINE RKMSG(IER,SRNAME,NREC,FLAG)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      To process error messages and terminate the program
C                in the event of a "catastrophic" failure.
C
C  Input:        IER, SRNAME, NREC
C  Output:       FLAG
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM7/ OUTCH
C                                /RKCOM8/ MSG, UTASK
C                                /RKCOM9/ REC
C                Alters:         none
C
C  Comments:
C  =========
C  The output variable FLAG is assigned the value of the input variable IER.
C
C  IER = -2  reports a successful call of the subroutine SRNAME and
C            indicates that special action is to be taken elsewhere
C            in the suite.  FLAG is set and a return is effected.
C
C  IER = 1   reports a successful call of the subroutine SRNAME.  FLAG
C            is set and a return is effected.
C
C  1 < IER < 911 and MSG = .TRUE.: a message of NREC records contained in
C            the array REC(*) is written to the standard output channel, 
C            OUTCH.  FLAG is set and a return is effected.
C
C  IER = 911 reports a "catastrophic" error was detected in SRNAME.  A
C            message is written to OUTCH regardless of the value of MSG and
C            normally the execution of the program is terminated.  The
C            execution is not terminated if the error is the result of an
C            indirect call to CT, RESET, or INTRP through UT (UTASK = .TRUE.).
C            Termination can be prevented by using the subroutine SOFTFL.
C
C  IER = 912 reports that a "catastrophic" error was detected earlier and
C            termination was prevented, but the user has failed to take
C            appropriate remedial action.  Execution is terminated.
C
C     .. Scalar Arguments ..
      INTEGER           FLAG, IER, NREC
      CHARACTER*(*)     SRNAME
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Common Block for Integrator Options ..
      LOGICAL           MSG, UTASK
      COMMON /RKCOM8/   MSG, UTASK
      SAVE   /RKCOM8/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      INTEGER           PLUS1
      LOGICAL           ASK, TELL
      PARAMETER         (PLUS1=1,ASK=.TRUE.,TELL=.FALSE.)
C     .. Local Scalars ..
      INTEGER           I
      LOGICAL           BADERR, OK, ON, UTCALL
C     .. External Subroutines ..
      EXTERNAL          CHKFL, RKSIT, SOFTFL
C     .. Executable Statements ..
C
C  Check where the call came from - if it is an indirect call from UT,
C  the run is not STOPped.
      UTCALL = (SRNAME.EQ.'RESET' .OR. SRNAME.EQ.'CT' .OR.
     &          SRNAME.EQ.'INTRP') .AND. UTASK
C
C  Check if can continue with integrator.
      OK = (SRNAME.EQ.'CT' .OR. SRNAME.EQ.'UT') .AND.
     &     (IER.EQ.2 .OR. IER.EQ.3 .OR. IER.EQ.4)
C
C  Check if program termination has been overridden.
      CALL SOFTFL(ASK,ON)
C
      IF ((MSG.AND.IER.GT.PLUS1) .OR. IER.GE.911) THEN
         WRITE (OUTCH,'(/A)') ' **'
         WRITE (OUTCH,'(A)') (REC(I),I=1,NREC)
         IF (IER.GE.911) THEN
            WRITE (OUTCH,'(A/A,A,A/A/)') 
     &' **',
     &' ** Catastrophic error detected in ', SRNAME, '.', 
     &' **'
            IF ((.NOT.(UTCALL.OR.ON).AND.IER.EQ.911) .OR.
     &          IER.EQ.912) THEN
               WRITE (OUTCH,'(A/A/A)') 
     &' **',
     &' ** Execution of your program is being terminated.',
     &' **'
               STOP
            END IF
         ELSE IF (OK) THEN
            WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)') 
     &' **',
     &' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.',
     &' ** You can continue integrating this problem.',
     &' **'
         ELSE
            WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)') 
     &' **',
     &' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.',
     &' ** You cannot continue integrating this problem.', 
     &' **'
         END IF
      END IF
      DO 20 I = NREC + 1, 10
         REC(I) = ' '
   20 CONTINUE
      FLAG = IER
C
C  TELL RKSIT the status of the routine associated with SRNAME
      CALL RKSIT(TELL,SRNAME,FLAG)
C
C  Indicate that a catastrophic error has been detected
      BADERR = FLAG .GE. 911
      CALL CHKFL(TELL,BADERR)
C
      RETURN
C
      END
      SUBROUTINE RKSIT(ASK,SRNAME,STATE)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      To save or enquire about the status of each
C                subprogram in the suite.
C
C  Input:        ASK, SRNAME
C  Input/output: STATE
C
C
C  Comments:
C  =========
C  SRNAME indicates which routine is of interest in the call to RKSIT.
C
C  If ASK=.FALSE., then the value of STATE (which as used is the error
C  flag) for the routine SRNAME is saved internally. This value of STATE
C  is usually positive.  There are two exceptions:
C     1. SRNAME='SETUP' and STATE=-1 indicates a completely new problem,
C        so all SAVEd states are cleared.
C     2. STATE=-2 is used by some routines in the suite to indicate
C        circumstances which require special action.
C
C  If ASK=.TRUE., then RKSIT first checks to see if there were any
C  catastrophic errors, that is, a SAVEd state has value 911. This should
C  happen only when the user has overridden program termination in the event
C  of catastrophic failures from routines in the package but has failed to
C  take appropriate action. If this is the case, then RKSIT returns a value
C  of STATE = 911 which forces a termination of execution inside RKMSG. If 
C  no catastrophic errors are flagged, then STATE returns the saved state 
C  value for the routine specified by SRNAME.
C
C     .. Scalar Arguments ..
      INTEGER           STATE
      LOGICAL           ASK
      CHARACTER*(*)     SRNAME
C     .. Parameters ..
      INTEGER           STATES, MINUS1
      PARAMETER         (STATES=7,MINUS1=-1)
C     .. Local Scalars ..
      INTEGER           I, NAME
C     .. Local Arrays ..
      INTEGER           SVSTA(STATES)
C     .. Save statement ..
      SAVE              SVSTA
C     .. Data statements ..
      DATA              SVSTA/STATES*MINUS1/
C     .. Executable Statements ..
C
      IF (SRNAME.EQ.'SETUP') THEN
         NAME = 1
      ELSE IF (SRNAME.EQ.'UT') THEN
         NAME = 2
      ELSE IF (SRNAME.EQ.'STAT') THEN
         NAME = 3
      ELSE IF (SRNAME.EQ.'GLBERR') THEN
         NAME = 4
      ELSE IF (SRNAME.EQ.'CT') THEN
         NAME = 5
      ELSE IF (SRNAME.EQ.'INTRP') THEN
         NAME = 6
      ELSE IF (SRNAME.EQ.'RESET') THEN
         NAME = 7
      ELSE
         NAME = 0
      END IF
C
C  (Re)initialize if SETUP is telling RKSIT to do so.
      IF (.NOT.ASK .AND. NAME.EQ.1 .AND. STATE.EQ.MINUS1) THEN
         DO 20 I = 1, STATES
            SVSTA(I) = MINUS1
   20    CONTINUE
         GO TO 60
      END IF
C
C  Check for 911 on exit from a previous call.
      IF (ASK) THEN
         DO 40 I = 1, STATES
            IF (SVSTA(I).EQ.911) THEN
               STATE = 911
               GO TO 60
            END IF
   40    CONTINUE
      END IF
C
      IF (ASK) THEN
         STATE = SVSTA(NAME)
      ELSE
         SVSTA(NAME) = STATE
      END IF
C
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE TRUERR(F,NEQ,Y,TOL,WEIGHT,ZY,ZYP,ZERROR,ZYNEW,ZERRES,
     &                  ZSTAGE,IER)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      Compute a running RMS measure of the true (global) error
C                for a general Runge-Kutta pair.
C
C
C  Input:        NEQ, Y(*), TOL, WEIGHT(*),
C  Input/output: ZY(*), ZYP(*), ZERROR(*)
C  Workspace:    ZYNEW(*), ZERRES(*), ZSTAGE(NEQ,*)
C  Output:       IER
C  External:     F
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM2/ T, HOLD
C                                /RKCOM5/ TOOSML, ORDER, NSEC
C                                /RKCOM7/ TINY
C                Alters:         /RKCOM6/ MAXERR, LOCMAX, GNFCN
C
C  Comments:
C  =========
C  A secondary integration is performed using a fraction of the step size 
C  of the primary integration. ZY(*) and ZYP(*) are the approximate solution
C  and first derivative of this secondary integration. ZERRES(*) contains the 
C  error estimates for the secondary integration. ZYNEW(*) and ZSTAGE(*,*) are
C  workspace for taking a step. The error assessment is computed using the
C  difference of the primary and secondary solutions at the primary
C  integration points as an estimate of the true error there.  The weights 
C  used are those of the error test of the primary integration. This error 
C  assessment is maintained in the vector ZERROR(*).  MAXERR and LOCMAX 
C  contain the maximum contribution to the assessment and its location,
C  respectively.  The number of calls to F is counted by GNFCN.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  TOL
      INTEGER           IER, NEQ
C     .. Array Arguments ..
      DOUBLE PRECISION  WEIGHT(*), Y(*), ZERRES(*), ZERROR(*),
     &                  ZSTAGE(NEQ,*), ZY(*), ZYNEW(*), ZYP(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Parameters ..
      DOUBLE PRECISION  PT1, TEN, DUMMY
      PARAMETER         (PT1=0.1D0,TEN=10.0D0,DUMMY=1.0D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  DIFF, ERRMAX, HMIN, HSEC, MXERLC, TSEC, ZLERR,
     &                  ZTEST1, ZTEST2
      INTEGER           ISTEP, L, LEVEL
      LOGICAL           LDUMMY, MAIN
C     .. Local Arrays ..
      DOUBLE PRECISION  DUMARR(1)
C     .. External Subroutines ..
      EXTERNAL          STEP
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, MAX
C     .. Executable Statements ..
      TSEC = T - HOLD
      HSEC = HOLD/DBLE(NSEC)
      HMIN = MAX(TINY,TOOSML*MAX(ABS(TSEC),ABS(T)))
      IF (ABS(HSEC).LT.HMIN) THEN
         IER = 6
         GO TO 120
      END IF
      ZTEST1 = TOL/DBLE(NSEC)
      ZTEST2 = TOL/TEN
      LEVEL = 0
C
C  The subroutine STEP is used to take a step.  In its use in the primary
C  integration provision is made for getting on scale in the first step.
C  In this situation only the subroutine might reduce the step size.  By
C  setting MAIN = .FALSE., the subroutine will take a step of the size input.
C  In this use of the subroutine, all items of the call list appearing after
C  MAIN are dummy variables.
C
C  Perform secondary integration.
      MAIN = .FALSE.
      LDUMMY = .FALSE.
      DO 60 ISTEP = 1, NSEC
C
C  Take a step.
         CALL STEP(F,NEQ,TSEC,ZY,ZYP,ZSTAGE,ZTEST1,HSEC,WEIGHT,ZYNEW,
     &             ZERRES,ZLERR,MAIN,DUMMY,DUMARR,LDUMMY)
C
C  The primary integration is using a step size of HUSED and the secondary
C  integration is using the smaller step size HSEC = HUSED/NSEC.  If steps
C  of this size were taken from the same starting point and the asymptotic
C  behavior were evident, the smaller step size would result in a local error
C  that is considerably smaller, namely by a factor of 1/(NSEC**(ORDER+1)).
C  If the two approximate solutions are close and TOLR is neither too large nor
C  too small, this should be approximately true.  The step size is chosen in
C  the primary integration so that the local error ERR is no larger than TOLR.
C  The local error, ZLERR, of the secondary integration is compared to TOLR in
C  an attempt to diagnose a secondary integration that is not rather more
C  accurate than the primary integration.
C
         IF (ZLERR.GE.ZTEST1) THEN
            LEVEL = 2
         ELSE IF (ZLERR.GT.ZTEST2) THEN
            LEVEL = LEVEL + 1
         END IF
         IF (LEVEL.GE.2) THEN
            IER = 6
            GO TO 120
         END IF
C
C  Advance TSEC and the dependent variables ZY(*) and ZYP(*).
         TSEC = T - DBLE(NSEC-ISTEP)*HSEC
         DO 20 L = 1, NEQ
            ZY(L) = ZYNEW(L)
   20    CONTINUE
C
         IF (FSAL) THEN
C
C  When FSAL = .TRUE., the derivative ZYP(*) is the last stage of the step.
            DO 40 L = 1, NEQ
               ZYP(L) = ZSTAGE(L,LSTSTG)
   40       CONTINUE
         ELSE
C
C  Call F to evaluate ZYP(*).
            CALL F(TSEC,ZY,ZYP)
            GNFCN = GNFCN + 1
         END IF
C
   60 CONTINUE
C
C  Update the maximum error seen, MAXERR, and its location, LOCMAX.
C  Use local variables ERRMAX and MXERLC.
C
      ERRMAX = MAXERR
      MXERLC = LOCMAX
      DO 80 L = 1, NEQ
         DIFF = ABS(ZY(L)-Y(L))/WEIGHT(L)
         IF (DIFF.GT.ERRMAX) THEN
            ERRMAX = DIFF
            MXERLC = T
         END IF
   80 CONTINUE
C
C  If the global error is greater than 0.1D0, the solutions have diverged so
C  far that comparing them may not provide a reliable estimate of the global
C  error. The test is made before ZERROR(*) and MAXERR, LCMXER are updated so
C  that on a failure, they refer to the last reliable results.
C
      IF (ERRMAX.GT.PT1) THEN
         IER = 6
         GO TO 120
      ELSE
         MAXERR = ERRMAX
         LOCMAX = MXERLC
         DO 100 L = 1, NEQ
            DIFF = ABS(ZY(L)-Y(L))/WEIGHT(L)
            ZERROR(L) = ZERROR(L) + DIFF**2
  100    CONTINUE
         IER = 1
      END IF
C
C  Exit point for TRUERR
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE STEP(F,NEQ,TNOW,Y,YP,STAGES,TOL,HTRY,WEIGHT,YNEW,
     &                ERREST,ERR,MAIN,HMIN,THRES,PHASE2)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      To compute a step of an explicit Runge-Kutta
C                method and estimate the local error of the step.
C
C  Input:        NEQ, TNOW, Y(*), YP(*), TOL, MAIN, HMIN, THRES(*)
C  Input/output: HTRY, PHASE2, LAST, WEIGHT(*)
C  Output:       STAGES(NEQ,*), YNEW(*), ERREST(*), ERR
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM1/ TND
C                                /RKCOM2/ LAST
C                                /RKCOM4/ A, B, C, BHAT, PTR, NSTAGE, METHD
C                                /RKCOM5/ FSAL
C                Alters:         /RKCOM2/ NFCN, LAST
C                                /RKCOM6/ GNFCN
C
C  Comments:
C  =========
C  From an approximate solution Y(*) at TNOW and first derivative there,
C  YP(*) = F(TNOW,Y,YP), a step is taken to get an approximation YNEW(*)
C  at TNOW + HTRY. The Runge-Kutta method and how it is used are defined
C  by A, B, C, BHAT, PTR, NSTAGE, METHD and FSAL. Intermediate stages
C  of the method are stored in the array STAGES(NEQ,*). The error in
C  each solution component is estimated and returned in ERREST(*). A
C  weighted maximum norm of the local error, ERR, is formed. For some
C  methods an intermediate error estimate can be computed before completion
C  of the step (see routine STEPB); if the estimate is greater than the
C  specified tolerance TOL, the computation of the step is terminated.
C
C  When global error estimation is desired, two integrations are done.
C  The usual integration is referred to as the "primary", or "main",
C  integration (MAIN=.TRUE.).  For global error estimation another,
C  "secondary" integration (MAIN=.FALSE.) is carried out with a smaller
C  step size.  The weight vector WEIGHT(*) used in computing ERR is
C  determined by the main integration.  Thus this argument is output when
C  MAIN = .TRUE. and input when MAIN = .FALSE..
C
C  When taking the first step in an integration, the logical variable
C  PHASE2 may be input as .TRUE. and if the first step is the whole of
C  the range of integration, then LAST will be .TRUE.. When PHASE2=.TRUE.,
C  the first three stages are monitored to help assure that the step
C  size H is small enough for the integration to be stable and for the
C  estimate of the error of the step to be credible. Calls are made to
C  the subroutine STEPA for this purpose. If necessary, H will be
C  reduced in STEPA (and LAST altered accordingly) and the step retried
C  in STEP until an acceptable value is found.
C
C  In the primary integration the number of calls to F is counted by
C  NFCN, and in the secondary integration, by GNFCN.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ERR, HMIN, HTRY, TNOW, TOL
      INTEGER           NEQ
      LOGICAL           MAIN, PHASE2
C     .. Array Arguments ..
      DOUBLE PRECISION  ERREST(*), STAGES(NEQ,*), THRES(*), WEIGHT(*),
     &                  Y(*), YNEW(*), YP(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Global Error Assessment ..
      DOUBLE PRECISION  MAXERR, LOCMAX
      INTEGER           GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
     &                  PRZYNU
      LOGICAL           ERASON, ERASFL
      COMMON /RKCOM6/   MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
     &                  PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
      SAVE   /RKCOM6/
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, HALF, ONE
      PARAMETER         (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  AVGY, TSTG
      INTEGER           I, J, L
      LOGICAL           CUTBAK
C     .. External Subroutines ..
      EXTERNAL          STEPA, STEPB
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SIGN
C     .. Executable Statements ..
C
C  Many of the following loops over L = 1, NEQ have constant array values
C  inside. The code is written with clarity in mind.  Any optimizing
C  compiler will identify these occurrences and take appropriate action.
C  A check for zero multipliers has been included so as to prevent
C  needless computation resulting from the storing of zero coefficients
C  in the arrays for the sake of clarity.  The array ERREST(*) is used
C  for working storage in this computation.
C
   20 CONTINUE
      IF (MAIN) THEN
         IF (PHASE2) THEN
C
C  Initialize weights for measuring the local error.
            DO 40 L = 1, NEQ
               WEIGHT(L) = MAX(THRES(L),ABS(Y(L)))
   40       CONTINUE
         END IF
      END IF
C
      DO 140 I = 2, NSTAGE
         DO 100 J = 1, I - 1
            IF (J.EQ.1) THEN
               DO 60 L = 1, NEQ
                  ERREST(L) = A(I,1)*YP(L)
   60          CONTINUE
            ELSE
               IF (A(I,J).NE.ZERO) THEN
                  DO 80 L = 1, NEQ
                     ERREST(L) = ERREST(L) + A(I,J)*STAGES(L,PTR(J))
   80             CONTINUE
               END IF
            END IF
  100    CONTINUE
         DO 120 L = 1, NEQ
            YNEW(L) = Y(L) + HTRY*ERREST(L)
  120    CONTINUE
C
C  METHD = 2 is special in that an estimate of the local error can be
C  formed before the step is completed.  If the step is a failure,
C  return immediately.  Otherwise, complete the step and compute a more
C  accurate error estimate.
         IF (METHD.EQ.2 .AND. I.EQ.7) THEN
            CALL STEPB(NEQ,Y,YP,HTRY,YNEW,STAGES,THRES,ERR,MAIN,WEIGHT)
            IF (ERR.GT.TOL) RETURN
         END IF
C
         TSTG = TNOW + C(I)*HTRY
         IF (MAIN .AND. LAST .AND. C(I).EQ.ONE) TSTG = TND
         CALL F(TSTG,YNEW,STAGES(1,PTR(I)))
C
C  Increment the counter for the number of function evaluations
C  depending on whether the primary or secondary integration is taking
C  place.
         IF (MAIN) THEN
            NFCN = NFCN + 1
         ELSE
            GNFCN = GNFCN + 1
         END IF
C
C----------------------------------------------------------------------
C  When PHASE2 is .TRUE. we are in the second phase of the automatic
C  selection of the initial step size.  The results of the first three
C  stages are monitored in the subroutine STEPA for evidence that H is
C  too large -- instability and/or an unreliable estimate of the error
C  of the step is then possible.  When the subroutine believes H to be
C  too large, it returns CUTBAK = .TRUE. and a suitably reduced H for
C  another try.
C
         IF (MAIN) THEN
            IF (PHASE2) THEN
               IF (I.LE.3 .AND. ABS(HTRY).GT.HMIN) THEN
                  CALL STEPA(TNOW,Y,YP,TSTG,YNEW,STAGES(1,PTR(I)),
     &                       HTRY,WEIGHT,CUTBAK)
                  IF (CUTBAK) THEN
                     LAST = .FALSE.
C
C  Make sure that STEPA does not reduce the step size below the
C  minimum. If it does, reset H to HMIN and deactivate PHASE2.
                     IF (ABS(HTRY).LE.HMIN) THEN
                        HTRY = SIGN(HMIN,HTRY)
                        PHASE2 = .FALSE.
                     END IF
                     GO TO 20
                  END IF
               END IF
            END IF
         END IF
C----------------------------------------------------------------------
C
  140 CONTINUE
C
C  Some formulas are constructed so that the last stage represents
C  the result of the step (FSAL=.TRUE.), hence if the step is acceptable,
C  it will be the first stage for the next step. When FSAL=.FALSE., we
C  have to complete the computation of the step.
C
      IF (.NOT.FSAL) THEN
         DO 200 I = 1, NSTAGE
            IF (I.EQ.1) THEN
               DO 160 L = 1, NEQ
                  ERREST(L) = BHAT(1)*YP(L)
  160          CONTINUE
            ELSE
               IF (BHAT(I).NE.ZERO) THEN
                  DO 180 L = 1, NEQ
                     ERREST(L) = ERREST(L) + BHAT(I)*STAGES(L,PTR(I))
  180             CONTINUE
               END IF
            END IF
  200    CONTINUE
         DO 220 L = 1, NEQ
            YNEW(L) = Y(L) + HTRY*ERREST(L)
  220    CONTINUE
      END IF
C
C  Form an estimate of the error in the lower order formula by comparing
C  it to the higher order formula of the pair. ERREST(*) has been used
C  as working storage above.  The higher order approximation has been
C  formed as YNEW(*) = Y(*) + HTRY*ERREST(*) where ERREST(*) is a linear
C  combination of the stages of the formula. The lower order result also
C  has the form Y(*) plus HTRY times a different linear combination of
C  the stages. Hence, this different linear combination of stages for
C  the lower order formula can just be subtracted from the combination
C  stored in ERREST(*) to produce the errors. The result is then
C  multiplied by HTRY to obtain the error estimate.
C
      DO 280 I = 1, NSTAGE
         IF (I.EQ.1 .AND. B(1).NE.ZERO) THEN
            DO 240 L = 1, NEQ
               ERREST(L) = ERREST(L) - B(1)*YP(L)
  240       CONTINUE
         ELSE
            IF (B(I).NE.ZERO) THEN
               DO 260 L = 1, NEQ
                  ERREST(L) = ERREST(L) - B(I)*STAGES(L,PTR(I))
  260          CONTINUE
            END IF
         END IF
  280 CONTINUE
      DO 300 L = 1, NEQ
         ERREST(L) = HTRY*ERREST(L)
  300 CONTINUE
C
C  The error in a solution component is measured relative to a weight
C  that is the larger of a threshold and the size of the solution over
C  the step.  Using the magnitude of a solution component at both ends
C  of the step in the definition of "size" increases the robustness of
C  the test. When global error estimation is specified, the weight
C  vector WEIGHT(*) is defined by the primary integration and is then
C  used in the secondary integration.
C
      IF (MAIN) THEN
         DO 320 L = 1, NEQ
            AVGY = HALF*(ABS(Y(L))+ABS(YNEW(L)))
            WEIGHT(L) = MAX(AVGY,THRES(L))
  320    CONTINUE
      END IF
C
      ERR = ZERO
      DO 340 L = 1, NEQ
         ERR = MAX(ERR,ABS(ERREST(L)/WEIGHT(L)))
  340 CONTINUE   
C
      RETURN
      END
      SUBROUTINE STEPA(TNOW,Y,YP,TSTG,YSTG,YPSTG,HTRY,WEIGHT,CUTBAK)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      To calculate an "on-scale" step size for phase 2 of
C                the initial step size computation.
C
C  Input:        TNOW, Y(*), YP(*), TSTG, YSTG(*), YPSTG(*)
C  Input/output: HTRY, WEIGHT
C  Output:       CUTBAK
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM1/ TND, NEQ
C                                /RKCOM5/ STBRAD, RS1, RS4
C                                /RKCOM7/ RNDOFF
C                Alters:         none
C
C  Comments:
C  =========
C  This subroutine is used during the first three stages of the first step.
C  A Lipschitz constant L for the differential equation in autonomous form
C  is approximated, and the product abs(HTRY)*L is compared to an approximate
C  radius, STBRAD, of the stability region of the method. The step size is 
C  reduced as necessary, within a range specified by the step size control 
C  parameters RS1 and RS4, to assure stability and give some confidence in 
C  the error estimator.  If HTRY is reduced, CUTBAK is set .TRUE..
C
C  Y(*) and YP(*) contain the solution and its derivative at TNOW and
C  similarly YSTG(*) and YPSTG(*) contain approximations at TSTG.
C
C  Normally the weights used in the control of the error depend on the
C  size of the solution at the beginning and at the end of the step, but
C  at this time we do not have a solution at the end of the step.  Each
C  stage YSTG(*) of the Runge - Kutta process represents a low order
C  approximation to the solution at TSTG.  Because the initial value of
C  WEIGHT(*) provided in the first phase of the scheme is based only on
C  the solution at T and THRES(*), it is continually updated in STEPA to
C  account for the size of the solution throughout the step as revealed
C  by the intermediate stages YSTG(*). Inside this subroutine only, the
C  differential equation is converted to autonomous form. After the
C  conversion, the end of the interval of integration, TND, is used
C  to define a suitable weight for the independent variable.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HTRY, TNOW, TSTG
      LOGICAL           CUTBAK
C     .. Array Arguments ..
      DOUBLE PRECISION  WEIGHT(*), Y(*), YP(*), YPSTG(*), YSTG(*)
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  ARGDIF, FDIFF, SCL, TDIFF, TWT, WT, YNRM, YSTGNM
      INTEGER           L
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MIN
C     .. Executable Statements ..
C
C  Update the weights to account for the current intermediate solution
C  approximation YSTG(*).  Compute the sizes of Y(*) and YSTG(*) in the
C  new norm.  The size of the Lipschitz constant is assessed by a difference
C  in the arguments Y(*), YSTG(*) and a difference in the function evaluated
C  at these arguments.
C
      YNRM = ZERO
      YSTGNM = ZERO
      ARGDIF = ZERO
      FDIFF = ZERO
      DO 20 L = 1, NEQN
         WT = MAX(WEIGHT(L),ABS(YSTG(L)))
         WEIGHT(L) = WT
         YNRM = MAX(YNRM,ABS(Y(L))/WT)
         YSTGNM = MAX(YSTGNM,ABS(YSTG(L))/WT)
         ARGDIF = MAX(ARGDIF,ABS(YSTG(L)-Y(L))/WT)
         FDIFF = MAX(FDIFF,ABS(YPSTG(L)-YP(L))/WT)
   20 CONTINUE
C
C  The transformation of the equation to autonomous form is done
C  implicitly.  The difference of the arguments must take into account
C  the difference between the values of the independent variable T and
C  TSTG. The difference of the corresponding component of the function
C  is zero because of the way the standard transformation is done.
C
      TDIFF = TSTG - TNOW
      TWT = ABS(TND-TNOW)
      YNRM = MAX(YNRM,ABS(TNOW)/TWT)
      YSTGNM = MAX(YSTGNM,ABS(TSTG)/TWT)
      ARGDIF = MAX(ARGDIF,ABS(TDIFF)/TWT)
C
C  The ratio FDIFF/ARGDIF is a lower bound for, and an approximation to, a
C  Lipschitz constant L for the differential equation written in autonomous
C  form.  First we must ask if the difference ARGDIF is significant in the 
C  precision available.  If it appears to be, we insist that abs(HTRY)*L be 
C  less than an approximate radius, STBRAD, of the stability region of the
C  method.  This is more stringent than necessary for stability, possibly a
C  lot more stringent, but the aim is to get an HTRY small enough that the
C  error estimate for the step is credible.  The reduction is required to be
C  at least as much as the step control parameter RS1. It is necessary to 
C  limit the reduction of HTRY at any one time because we may be misled in 
C  the size of the reduction that is appropriate due to nonlinearity of the 
C  differential equation and to inaccurate weights caused by HTRY much too 
C  large.  The reduction is not permitted to be more than the step control 
C  parameter RS4.
C
      CUTBAK = .FALSE.
      IF (ARGDIF.GT.RNDOFF*MAX(YNRM,YSTGNM)) THEN
         IF ((ABS(HTRY)*FDIFF).GT.(STBRAD*ARGDIF)) THEN
            SCL = (STBRAD*ARGDIF)/(ABS(HTRY)*FDIFF)
            SCL = MIN(SCL,RS1)
            SCL = MAX(SCL,RS4)
            HTRY = SCL*HTRY
            CUTBAK = .TRUE.
         END IF
      END IF
C
      RETURN
      END
      SUBROUTINE STEPB(NEQ,Y,YP,H,YNEW,STAGES,THRES,ERR,MAIN,WEIGHT)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      To compute an error estimate for METHD = 2 prior
C                to completing the step.
C
C  Input:        NEQ, Y(*), YP(*), H, STAGES(NEQ,*), THRES(*), MAIN,
C                WEIGHT(*)
C  Output:       ERR
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM4/ E, PTR
C                Alters:         none
C
C  Comments:
C  =========
C  If global error assessment is taking place, then MAIN = .FALSE. and
C  the weight vector generated by the primary integration is used.  The
C  error estimate is a linear combination (with coefficients in E(*))
C  of the stages stored in STAGES(*,*) (located by PTR(*)).
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ERR, H
      INTEGER           NEQ
      LOGICAL           MAIN
C     .. Array Arguments ..
      DOUBLE PRECISION  STAGES(NEQ,*), THRES(*), WEIGHT(*), Y(*),
     &                  YNEW(*), YP(*)
C     .. Common Block to hold Formula Definitions ..
      DOUBLE PRECISION  A(13,13), B(13), C(13), BHAT(13), R(11,6),
     &                  E(7)
      INTEGER           PTR(13), NSTAGE, METHD, MINTP
      LOGICAL           INTP
      COMMON /RKCOM4/   A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
     &                  MINTP, INTP
      SAVE   /RKCOM4/
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, HALF
      PARAMETER         (ZERO=0.0D0,HALF=0.5D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  AVGY, SUM, WT
      INTEGER           INDEX, L
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX
C     .. Executable Statements ..
C
      ERR = ZERO
      DO 40 L = 1, NEQ
C
C  Estimate the local error of component L. The coding makes use of
C  E(2) = 0.0D0 and E(7) = 0.0D0.
C
         SUM = E(1)*YP(L)
         DO 20 INDEX = 3, 6
            SUM = SUM + E(INDEX)*STAGES(L,PTR(INDEX))
   20    CONTINUE
C
C  The local error is H*SUM.  A weighted maximum norm of SUM is formed 
C  and then the factor of H is taken into account.
C  
         IF (MAIN) THEN
            AVGY = HALF*(ABS(Y(L))+ABS(YNEW(L)))
            WT = MAX(AVGY,THRES(L))
         ELSE
            WT = WEIGHT(L)
         END IF
C
         ERR = MAX(ERR,ABS(SUM/WT))
   40 CONTINUE
      ERR = ABS(H)*ERR
C
      RETURN
      END
      SUBROUTINE STIFF(F,HAVG,JFLSTP,TOOMCH,MAXFCN,WORK,IER,NREC)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:      Diagnose stiffness.  This depends on two things: whether
C                the step size is being restricted on grounds of stability
C                and whether the integration to TND can be completed in no
C                more than MAXFCN function evaluations.
C
C  Input:        HAVG, TOOMCH, MAXFCN, WORK(*)
C  Input/output: JFLSTP
C  Output:       IER, NREC
C  Workspace:    WORK(*)
C  External:     F
C
C  Common:       Initializes:    /RKCOM9/ REC
C                Reads:          /RKCOM1/ TND, NEQN
C                                /RKCOM2/ T, H, NFCN, SVNFCN, OKSTP
C                                /RKCOM3/ PRY, PRYP, PRTHRS, PRWT, PRSCR,
C                                         PRSTGS, PRYOLD
C                                /RKCOM5/ COST
C                Alters:         /RKCOM2/ NFCN
C                                /RKCOM9/ REC
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HAVG
      INTEGER           IER, JFLSTP, MAXFCN, NREC
      LOGICAL           TOOMCH
C     .. Array Arguments ..
      DOUBLE PRECISION  WORK(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Common Block for General Workspace Pointers ..
      INTEGER           PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      COMMON /RKCOM3/   PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
     &                  PRSTGS, PRINTP, LNINTP
      SAVE   /RKCOM3/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Error Message ..
      CHARACTER*80      REC(10)
      COMMON /RKCOM9/   REC
      SAVE   /RKCOM9/
C     .. Parameters ..
      DOUBLE PRECISION  HALF
      PARAMETER         (HALF=0.5D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  AVGY, XTRAWK
      INTEGER           L
      LOGICAL           LOTSFL, STIF, UNSURE
C     .. External Subroutines ..
      EXTERNAL          STIFFA
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, MAX, MOD
C     .. Executable Statements ..
C
      IF (MOD(OKSTP-10,40).EQ.0) THEN
         LOTSFL = JFLSTP .GE. 10
         JFLSTP = 0
      ELSE
         LOTSFL = .FALSE.
      END IF
C
C  If either too much work has been done or there are lots of failed steps,
C  test for stiffness.
C
      IF (TOOMCH .OR. LOTSFL) THEN
C
C  Regenerate weight vector
         DO 20 L = 1, NEQN
            AVGY = HALF*(ABS(WORK(PRY-1+L))+ABS(WORK(PRYOLD-1+L)))
            WORK(PRWT-1+L) = MAX(AVGY,WORK(PRTHRS-1+L))
   20    CONTINUE
C
C  STIFFA determines whether the problem is STIFF. In some circumstances it
C  is UNSURE.  The decision depends on two things: whether the step size is
C  being restricted on grounds of stability and whether the integration to
C  TND can be completed in no more than MAXFCN function evaluations.  The
C  last four arguments of STIFFA are vectors of length NEQN used for working
C  storage.  Some storage in WORK(*) reserved for the stages (there are a
C  minimum of three such vectors reserved for the METHDs implemented) and
C  the scratch vector starting at PRSCR are used for this purpose.
C
         CALL STIFFA(F,T,WORK(PRY),H,HAVG,TND,MAXFCN,WORK(PRWT),
     &               WORK(PRYP),WORK(PRERST),UNSURE,STIF,WORK(PRSTGS),
     &               WORK(PRSTGS+NEQN),WORK(PRSTGS+2*NEQN),WORK(PRSCR))
         IF (.NOT.UNSURE) THEN
            IF (STIF) THEN
C
C  Predict how much eXTRA WorK will be needed to reach TND.
               XTRAWK = (COST*ABS((TND-T)/HAVG))/DBLE(SVNFCN+NFCN)
               IER = 4
               WRITE (REC(NREC+1),'(A)')
     &' ** Your problem has been diagnosed as stiff.  If the '
               WRITE (REC(NREC+2),'(A,D13.5)')
     &' ** situation persists, it will cost roughly ', XTRAWK
               WRITE (REC(NREC+3),'(A)')
     &' ** times as much to reach TEND as it has cost to reach TNOW.'
               WRITE (REC(NREC+4),'(A)')
     &' ** You should probably change to a code intended for '
               WRITE (REC(NREC+5),'(A)')
     &' ** stiff problems. '
               NREC = NREC + 5
            END IF
         END IF
      END IF
C
      RETURN
      END
      SUBROUTINE STIFFA(F,X,Y,HNOW,HAVG,XEND,MAXFCN,WT,FXY,V0,UNSURE,
     &                  STIF,V1,V2,V3,VTEMP)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  External:     F
C  Input:        X, Y(*), HNOW, HAVG, XEND, MAXFCN, WT(*), FXY(*)
C  Input/Output  V0(*)
C  Output:       UNSURE, STIF
C  Workspace:    V1(*), V2(*), V3(*), VTEMP(*)
C
C  Common:       Initializes:    none
C                Reads:          /RKCOM1/ TND, NEQN
C                                /RKCOM5/ COST, STBRAD, TANANG
C                                /RKCOM7/ SQRRMC, CUBRMC
C                Alters:         none
C
C  STIFFA diagnoses stiffness for an explicit Runge-Kutta code.  When it
C  is called, either many step failures have been observed, or a lot of
C  work has been done.
C
C  The NEQ equations of the problem are defined by the subroutine F(X,Y,YP).
C  When STIFFA is called, the integration has reached X where the approximate
C  solution is Y(*).  The vector FXY(*) is defined by a call of F(X,Y,FXY).
C  It is an input argument because it is usually available from the integrator.
C
C  The last successful step was of size HNOW, and an average step size is
C  HAVG.  A weighted norm is used to measure the local error with the error
C  in solution component L divided by the positive weight WT(L) provided in 
C  the vector WT(*).
C
C  Explicit Runge - Kutta codes estimate the local error of Y(*) by
C  forming the difference of two approximate solutions.  This difference
C  must be provided in the vector V0(*).  When this difference is too
C  small to be significant, STIFFA will replace it with a "random" vector.
C
C  STIF is set .TRUE. when the average step size appears to be restricted
C  on grounds of stability.  In certain cases the variable UNSURE is set 
C  .TRUE.; the value of STIF is then not defined.
C
C  The stability region of the explicit Runge-Kutta formula is described
C  by quantities TANANG and STBRAD that are communicated by the setup routine
C  via COMMON.  Stability regions often change sharply near the imaginary
C  axis so that it is difficult to classify the stiffness of a problem with
C  eigenvalues of a local Jacobian that are "near" the imaginary axis.  For
C  this reason,  we consider only points Z in the upper left half complex
C  plane for which TAN( IMAG(Z)/( - RE(Z))) <= TANANG. Eigenvalues outside
C  this region are one reason for the code being UNSURE.  The stability
C  region is approximated by the intersection of a disk with this sector.
C  The radius of this disk is called STBRAD.
C
C  Working storage must be provided via the four vectors V1(*),V2(*),
C  V3(*),VTEMP(*).  These vectors must be of length at least NEQ.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HAVG, HNOW, X, XEND
      INTEGER           MAXFCN
      LOGICAL           STIF, UNSURE
C     .. Array Arguments ..
      DOUBLE PRECISION  FXY(*), V0(*), V1(*), V2(*), V3(*), VTEMP(*),
     &                  WT(*), Y(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Formula Characterisitcs ..
      DOUBLE PRECISION  TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4
      INTEGER           ORDER, LSTSTG, MAXTRY, NSEC
      LOGICAL           FSAL
      COMMON /RKCOM5/   TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
     &                  RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
     &                  NSEC, FSAL
      SAVE   /RKCOM5/
C     .. Common Block for Environment Parameters ..
      DOUBLE PRECISION  MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
      INTEGER           OUTCH
      COMMON /RKCOM7/   MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
     &                  OUTCH
      SAVE   /RKCOM7/
C     .. Parameters ..
      DOUBLE PRECISION  LARGE
      PARAMETER         (LARGE=1.0D+10)
      DOUBLE PRECISION  ZERO, P001, P9, ONE, TWO, FIVE, FIFTH
      PARAMETER         (ZERO=0.0D+0,P001=0.001D+0,P9=0.9D+0,ONE=1.0D+0,
     &                  TWO=2.0D+0,FIVE=5.0D+0,FIFTH=0.2D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  ALPHA1, ALPHA2, BETA1, BETA2, D1, D2, DET1,
     &                  DET2, DIST, RES2, RHO, RHO2, ROLD, SCALE, V0NRM,
     &                  V0V0, V0V1, V0V2, V1V1, V1V2, V1V3, V2V2, V2V3,
     &                  V3NRM, V3V3, XTRFCN, YNRM
      INTEGER           L, NTRY
      LOGICAL           ROOTRE
C     .. Local Arrays ..
      DOUBLE PRECISION  R1(2), R2(2), ROOT1(2), ROOT2(2)
C     .. External Functions ..
      DOUBLE PRECISION  DOTPRD
      EXTERNAL          DOTPRD
C     .. External Subroutines ..
      EXTERNAL          STIFFB, STIFFC, STIFFD
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MIN, SQRT
C     .. Executable Statements ..
C
C  If the current step size differs substantially from the average,
C  the problem is not stiff.
C
      IF (ABS(HNOW/HAVG).GT.FIVE .OR. ABS(HNOW/HAVG).LT.FIFTH) THEN
         STIF = .FALSE.
         UNSURE = .FALSE.
         RETURN
      ELSE
         UNSURE = .TRUE.
      END IF
C
C  The average step size is used to predict the cost in function evaluations
C  of finishing the integration to XEND.  If this cost is no more than MAXFCN,
C  the problem is declared not stiff: If the step size is being restricted on
C  grounds of stability, it will stay close to HAVG.  The prediction will
C  then be good, but the cost is too low to consider the problem stiff.  If
C  the step size is not close to HAVG, the problem is not stiff.  Either way
C  there is no point to testing for a step size restriction due to stability.
C
      XTRFCN = COST*ABS((XEND-X)/HAVG)
      IF (XTRFCN.LE.MAXFCN) THEN
         STIF = .FALSE.
         UNSURE = .FALSE.
         RETURN
      ELSE
         UNSURE = .TRUE.
      END IF
C
C  There have been many step failures or a lot of work has been done.  Now 
C  we must determine if this is due to the stability characteristics of the
C  formula.  This is done by calculating the dominant eigenvalues of the
C  local Jacobian and then testing whether HAVG corresponds to being on the
C  boundary of the stability region.
C
C  The size of Y(*) provides scale information needed to approximate
C  the Jacobian by differences.
C
      YNRM = SQRT(DOTPRD(Y,Y,WT,NEQN))
      SCALE = YNRM*SQRRMC
      IF (SCALE.EQ.ZERO) THEN
C
C  Degenerate case.  Y(*) is (almost) the zero vector so the scale is not 
C  defined.  The input vector V0(*) is the difference between Y(*) and a 
C  lower order approximation to the solution that is within the error 
C  tolerance.  When Y(*) vanishes, V0(*) is itself an acceptable approximate
C  solution, so we take SCALE from it, if this is possible.
C
         YNRM = SQRT(DOTPRD(V0,V0,WT,NEQN))
         SCALE = YNRM*SQRRMC
         IF (SCALE.EQ.ZERO) THEN
            UNSURE = .TRUE.
            RETURN
         END IF
      END IF
C
      V0V0 = DOTPRD(V0,V0,WT,NEQN)
      IF (V0V0.EQ.ZERO) THEN
C
C  Degenerate case.  V0(*) is (almost) the zero vector so cannot
C  be used to define a direction for an increment to Y(*).  Try a
C  "random" direction.
C
         DO 20 L = 1, NEQN
            V0(L) = ONE
   20    CONTINUE
         V0V0 = DOTPRD(V0,V0,WT,NEQN)
      END IF
      V0NRM = SQRT(V0V0)
      DO 40 L = 1, NEQN
         V0(L) = V0(L)/V0NRM
   40 CONTINUE
      V0V0 = ONE
C
C  Use a nonlinear power method to estimate the two dominant eigenvalues.
C  V0(*) is often very rich in the two associated eigenvectors.  For this 
C  reason the computation is organized with the expectation that a minimal 
C  number of iterations will suffice.  Indeed, it is necessary to recognize 
C  a kind of degeneracy when there is a dominant real eigenvalue.  The
C  subroutine STIFFB does this.  In the first try, NTRY = 1, a Rayleigh 
C  quotient for such an eigenvalue is initialized as ROLD.  After each 
C  iteration, REROOT computes a new Rayleigh quotient and tests whether the
C  two approximations agree to one tenth of one per cent and the eigenvalue,
C  eigenvector pair satisfy a stringent test on the residual.  ROOTRE = .TRUE.
C  signals that a single dominant real root has been found.
C
      NTRY = 1
   60 CONTINUE
C
      CALL STIFFD(V0,HAVG,X,Y,F,FXY,WT,SCALE,V0V0,V1,V1V1,VTEMP)
C
C  The quantity SQRT(V1V1/V0V0) is a lower bound for the product of HAVG
C  and a Lipschitz constant.  If it should be LARGE, stiffness is not
C  restricting the step size to the stability region.  The principle is
C  clear enough, but the real reason for this test is to recognize an
C  extremely inaccurate computation of V1V1 due to finite precision
C  arithmetic in certain degenerate circumstances.
C
      IF (SQRT(V1V1).GT.LARGE*SQRT(V0V0)) THEN
         UNSURE = .TRUE.
         RETURN
      END IF
C
      V0V1 = DOTPRD(V0,V1,WT,NEQN)
      IF (NTRY.EQ.1) THEN
         ROLD = V0V1/V0V0
C
C  This is the first Rayleigh quotient approximating the product of HAVG
C  and a dominant real eigenvalue.  If it should be very small, the
C  problem is not stiff.  It is important to test for this possibility so
C  as to prevent underflow and degeneracies in the subsequent iteration.
C
         IF (ABS(ROLD).LT.CUBRMC) THEN
            UNSURE = .FALSE.
            STIF = .FALSE.
            RETURN
         END IF
      ELSE
         CALL STIFFB(V1V1,V0V1,V0V0,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
         IF (ROOTRE) GO TO 100
      END IF
      CALL STIFFD(V1,HAVG,X,Y,F,FXY,WT,SCALE,V1V1,V2,V2V2,VTEMP)
      V0V2 = DOTPRD(V0,V2,WT,NEQN)
      V1V2 = DOTPRD(V1,V2,WT,NEQN)
      CALL STIFFB(V2V2,V1V2,V1V1,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
      IF (ROOTRE) GO TO 100
C
C  Fit a quadratic in the eigenvalue to the three successive iterates
C  V0(*),V1(*),V2(*) of the power method to get a first approximation to
C  a pair of eigenvalues.  A test made earlier in STIFFB implies that
C  the quantity DET1 here will not be too small.
C
      DET1 = V0V0*V1V1 - V0V1**2
      ALPHA1 = (-V0V0*V1V2+V0V1*V0V2)/DET1
      BETA1 = (V0V1*V1V2-V1V1*V0V2)/DET1
C
C  Iterate again to get V3, test again for degeneracy, and then fit a
C  quadratic to V1(*),V2(*),V3(*) to get a second approximation to a pair
C  of eigenvalues.
C
      CALL STIFFD(V2,HAVG,X,Y,F,FXY,WT,SCALE,V2V2,V3,V3V3,VTEMP)
      V1V3 = DOTPRD(V1,V3,WT,NEQN)
      V2V3 = DOTPRD(V2,V3,WT,NEQN)
      CALL STIFFB(V3V3,V2V3,V2V2,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
      IF (ROOTRE) GO TO 100
      DET2 = V1V1*V2V2 - V1V2**2
      ALPHA2 = (-V1V1*V2V3+V1V2*V1V3)/DET2
      BETA2 = (V1V2*V2V3-V2V2*V1V3)/DET2
C
C  First test the residual of the quadratic fit to see if we might
C  have determined a pair of eigenvalues.
C
      RES2 = ABS(V3V3+V2V2*ALPHA2**2+V1V1*BETA2**2+TWO*V2V3*ALPHA2+
     &       TWO*V1V3*BETA2+TWO*V1V2*ALPHA2*BETA2)
      IF (RES2.LE.V3V3*P001**2) THEN
C
C  Calculate the two approximate pairs of eigenvalues.
C
         CALL STIFFC(ALPHA1,BETA1,R1,R2)
         CALL STIFFC(ALPHA2,BETA2,ROOT1,ROOT2)
C
C  The test for convergence is done on the larger root of the second
C  approximation.  It is complicated by the fact that one pair of roots 
C  might be real and the other complex.  First calculate the spectral 
C  radius RHO of HAVG*J as the magnitude of ROOT1.  Then see if one of 
C  the roots R1,R2 is within one per cent of ROOT1.  A subdominant root 
C  may be very poorly approximated if its magnitude is much smaller than 
C  RHO -- this does not matter in our use of these eigenvalues.
C
         RHO = SQRT(ROOT1(1)**2+ROOT1(2)**2)
         D1 = (ROOT1(1)-R1(1))**2 + (ROOT1(2)-R1(2))**2
         D2 = (ROOT1(1)-R2(1))**2 + (ROOT1(2)-R2(2))**2
         DIST = SQRT(MIN(D1,D2))
         IF (DIST.LE.P001*RHO) GO TO 100
      END IF
C
C  Do not have convergence yet.  Because the iterations are cheap, and
C  because the convergence criterion is stringent, we are willing to try
C  a few iterations.
C
      IF (NTRY.LT.MAXTRY) THEN
         NTRY = NTRY + 1
         V3NRM = SQRT(V3V3)
         DO 80 L = 1, NEQN
            V0(L) = V3(L)/V3NRM
   80    CONTINUE
         V0V0 = ONE
         GO TO 60
      ELSE
         UNSURE = .TRUE.
         RETURN
      END IF
C
C                        **************
C
C  We now have the dominant eigenvalues.  Decide if the average step
C  size is being restricted on grounds of stability.  Check the real
C  parts of the eigenvalues.  First see if the dominant eigenvalue is
C  in the left half plane -- there won't be a stability restriction
C  unless it is. If there is another eigenvalue of comparable magnitude
C  with a positive real part, the problem is not stiff. If the dominant
C  eigenvalue is too close to the imaginary axis, we cannot diagnose
C  stiffness.
C
  100 CONTINUE
      IF (ROOT1(1).GT.ZERO) THEN
         STIF = .FALSE.
         UNSURE = .FALSE.
         RETURN
      END IF
      RHO2 = SQRT(ROOT2(1)**2+ROOT2(2)**2)
      IF (RHO2.GE.P9*RHO .AND. ROOT2(1).GT.ZERO) THEN
         STIF = .FALSE.
         UNSURE = .FALSE.
         RETURN
      END IF
      IF (ABS(ROOT1(2)).GT.ABS(ROOT1(1))*TANANG) THEN
         UNSURE = .TRUE.
         RETURN
      END IF
C
C  If the average step size corresponds to being well within the
C  stability region, the step size is not being restricted because
C  of stability.
C
      STIF = RHO .GE. P9*STBRAD
      UNSURE = .FALSE.
      RETURN
      END
      SUBROUTINE STIFFB(V1V1,V0V1,V0V0,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Input:        V1V1, V0V1, V0V0
C  Input/output: ROLD
C  Output:       RHO, ROOT1(*),ROOT2(*),ROOTRE
C
C  Decide if the iteration has degenerated because of a strongly
C  dominant real eigenvalue.  Have just computed the latest iterate.
C  V1V1 is its dot product with itself, V0V1 is the dot product
C  of the previous iterate with the current one, and V0V0 is the
C  dot product of the previous iterate with itself.  ROLD is a
C  previous Rayleigh quotient approximating a dominant real
C  eigenvalue.  It must be computed directly the first time the
C  subroutine is called.  It is updated each call to STIFFB, hence
C  is available for subsequent calls.
C
C  If there is a strongly dominant real eigenvalue, ROOTRE is set
C  .TRUE., ROOT1(*) returns the eigenvalue, RHO returns the magnitude
C  of the eigenvalue, and ROOT2(*) is set to zero.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  RHO, ROLD, V0V0, V0V1, V1V1
      LOGICAL           ROOTRE
C     .. Array Arguments ..
      DOUBLE PRECISION  ROOT1(2), ROOT2(2)
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, P001
      PARAMETER         (ZERO=0.0D+0,P001=0.001D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  DET, R, RES
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Executable Statements ..
C
      R = V0V1/V0V0
      RHO = ABS(R)
      DET = V0V0*V1V1 - V0V1**2
      RES = ABS(DET/V0V0)
      ROOTRE = DET .EQ. ZERO .OR. (RES.LE.V1V1*P001**2 .AND.
     &         ABS(R-ROLD).LE.P001*RHO)
      IF (ROOTRE) THEN
         ROOT1(1) = R
         ROOT1(2) = ZERO
         ROOT2(1) = ZERO
         ROOT2(2) = ZERO
      END IF
      ROLD = R
C
      RETURN
      END
      SUBROUTINE STIFFC(ALPHA,BETA,R1,R2)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Input:  ALPHA, BETA
C  Output: R1(*), R2(*)
C
C  This subroutine computes the two complex roots R1 and R2 of
C  the quadratic equation X**2 + ALPHA*X + BETA = 0.  The magnitude
C  of R1 is greater than or equal to the magnitude of R2. R1 and R2 are
C  returned as vectors of two components with the first being the real
C  part of the complex number and the second being the imaginary part.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ALPHA, BETA
C     .. Array Arguments ..
      DOUBLE PRECISION  R1(2), R2(2)
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, TWO
      PARAMETER         (ZERO=0.0D+0,TWO=2.0D+0)
C     .. Local Scalars ..
      DOUBLE PRECISION  DISC, SQDISC, TEMP
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SQRT
C     .. Executable Statements ..
      TEMP = ALPHA/TWO
      DISC = TEMP**2 - BETA
      IF (DISC.EQ.ZERO) THEN
C
C  Double root.
C
         R1(1) = -TEMP
         R1(2) = ZERO
         R2(1) = R1(1)
         R2(2) = R1(2)
         RETURN
      END IF
C
      SQDISC = SQRT(ABS(DISC))
      IF (DISC.LT.ZERO) THEN
C
C  Complex conjugate roots.
C
         R1(1) = -TEMP
         R1(2) = SQDISC
         R2(1) = R1(1)
         R2(2) = -R1(2)
      ELSE
C
C  Real pair of roots.  Calculate the bigger one in R1(1).
C
         IF (TEMP.GT.ZERO) THEN
            R1(1) = -TEMP - SQDISC
         ELSE
            R1(1) = -TEMP + SQDISC
         END IF
         R1(2) = ZERO
         R2(1) = BETA/R1(1)
         R2(2) = ZERO
      END IF
C
      RETURN
      END
      SUBROUTINE STIFFD(V,HAVG,X,Y,F,FXY,WT,SCALE,VDOTV,Z,ZDOTZ,VTEMP)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  External:     F
C  Input:        V(*), HAVG, X, Y(*), FXY(*), WT(*), SCALE, VDOTV,
C  Output:       Z(*), ZDOTZ
C  Workspace:    VTEMP(*)
C
C  For an input vector V(*) of length NEQ, this subroutine computes a vector
C  Z(*) that approximates the product HAVG*J*V where HAVG is an input scalar 
C  and J is the Jacobian matrix of a function F evaluated at the input 
C  arguments (X,Y(*)).  This function is defined by a subroutine of the form
C  F(T,U,F) that when given T and U(*), returns the value of the function in 
C  F(*).  The input vector FXY(*) is defined by F(X,Y,FXY).  Scaling is a 
C  delicate matter.  A weighted Euclidean norm is used with the (positive) 
C  weights provided in WT(*).  The input scalar SCALE is the square root of 
C  the unit roundoff times the norm of Y(*).  The square of the norm of the
C  input vector V(*) is input as VDOTV.  The routine outputs the square of
C  the norm of the output vector Z(*) as ZDOTZ.  The subroutine calls the
C  DOUBLE PRECISION FUNCTION DOTPRD(U,V,WT,NEQ) to compute the dot (inner)
C  product.  The vector VTEMP(*) is used for working storage.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HAVG, SCALE, VDOTV, X, ZDOTZ
C     .. Array Arguments ..
      DOUBLE PRECISION  FXY(*), V(*), VTEMP(*), WT(*), Y(*), Z(*)
C     .. Subroutine Arguments ..
      EXTERNAL          F
C     .. Common Block for Problem Definition ..
      DOUBLE PRECISION  TSTRT, TND, DIR, HSTRT, TOLR
      INTEGER           NEQN
      COMMON /RKCOM1/   TSTRT, TND, DIR, HSTRT, TOLR, NEQN
      SAVE   /RKCOM1/
C     .. Common Block to hold Problem Status ..
      DOUBLE PRECISION  T, H, TOLD, HOLD
      INTEGER           NFCN, SVNFCN, OKSTP, FLSTP
      LOGICAL           FIRST, LAST
      COMMON /RKCOM2/   T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
     &                  FIRST, LAST
      SAVE   /RKCOM2/
C     .. Local Scalars ..
      DOUBLE PRECISION  TEMP1, TEMP2
      INTEGER           L
C     .. External Functions ..
      DOUBLE PRECISION  DOTPRD
      EXTERNAL          DOTPRD
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
C
C  Scale V(*) so that it can be used as an increment to Y(*)
C  for an accurate difference approximation to the Jacobian.
C
      TEMP1 = SCALE/SQRT(VDOTV)
      DO 20 L = 1, NEQN
         VTEMP(L) = Y(L) + TEMP1*V(L)
   20 CONTINUE
C
      CALL F(X,VTEMP,Z)
      NFCN = NFCN + 1
C
C  Form the difference approximation.  At the same time undo
C  the scaling of V(*) and introduce the factor of HAVG.
C
      TEMP2 = HAVG/TEMP1
      DO 40 L = 1, NEQN
         Z(L) = TEMP2*(Z(L)-FXY(L))
   40 CONTINUE
C
      ZDOTZ = DOTPRD(Z,Z,WT,NEQN)
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION DOTPRD(U,V,WT,NEQ)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C  Purpose:   To compute a weighted Euclidean dot (inner) product of
C             two vectors.
C
C  Input:     U(*), V(*), WT(*), NEQ
C  Output:    the result DOTPRD is returned via the subprogram name
C
C  Comments:
C  =========
C  The vectors U(*), V(*), and WT(*) are of length NEQ. The components
C  of WT(*) are weights that must be non-zero.
C
C     .. Scalar Arguments ..
      INTEGER           NEQ
C     .. Array Arguments ..
      DOUBLE PRECISION  U(*), V(*), WT(*)
C     .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D0)
C     .. Local Scalars ..
      DOUBLE PRECISION  SUM
      INTEGER           L
C     .. Executable Statements ..
C
      SUM = ZERO
      DO 20 L = 1, NEQ
         SUM = SUM + (U(L)/WT(L))*(V(L)/WT(L))
   20 CONTINUE
C
      DOTPRD = SUM
C
      RETURN
      END
      SUBROUTINE SOFTFL(ASK,ON)
C
C  Purpose:      To prevent a program STOP after a "catastrophic"
C                failure when using a routine from RKSUITE.
C
C  Input:        ASK
C  Input/output: ON
C
C  Comments:
C  =========
C  When a "catastrophic" failure is detected, the default action of
C  RKSUITE is to write an explanation to the standard output channel,
C  OUTCH, and STOP.  This subroutine can be used to prevent the STOP and
C  so allow the main program to continue.  To do this, you call SOFTFL with
C  ASK = .FALSE. and ON = .TRUE.  You must then call the subroutine CHKFL
C  after every call to a user-callable routine in RKSUITE to check whether
C  a catastrophic error occurred and take appropriate action if it did.  Of
C  course, you may call SETUP at any time to start a new problem, but calling
C  any other user-callable routine in RKSUITE after a catastrophic error will
C  lead to a STOP (even when "soft failure" has been set "on").
C
C  When ON is set by a call to SOFTFL with ASK = .FALSE., the value of ON
C  is SAVEd.  The subroutine RKMSG in RKSUITE calls SOFTFL with ASK = .TRUE.
C  to find out the SAVEd value of ON.
C
C     .. Scalar Arguments ..
      LOGICAL           ASK, ON
C     .. Local Scalars ..
      LOGICAL           SOFT
C     .. Save statement ..
      SAVE              SOFT
C     .. Data statements ..
      DATA              SOFT/.FALSE./
C     .. Executable Statements ..
C
      IF (ASK) THEN
         ON = SOFT
      ELSE
         SOFT = ON
      END IF
C
      RETURN
      END
      SUBROUTINE CHKFL(ASK,ERROR)
C
C  Purpose:      Enquiry routine used in conjunction with SOFTFL.
C                Reports whether a "catastrophic" error was detected.
C
C  Input:        ASK
C  Input/output: ERROR
C
C  Comments:
C  =========
C  When a "catastrophic" failure is detected, the default action of
C  RKSUITE is to write an explanation to the standard output channel,
C  OUTCH, and STOP.  SOFTFL can be used to prevent the STOP and so
C  allow the main program to continue.  It is then necessary to call
C  CHKFL with ASK = .TRUE. after every call to a user-callable routine 
C  in RKSUITE to check whether a catastrophic error occurred and take 
C  appropriate action if it did.  If there was a catastrophic error, 
C  ERROR is returned .TRUE.  Of course, you may call SETUP at any time 
C  to start a new problem, but calling any other user-callable routine 
C  in RKSUITE after a catastrophic error will lead to a STOP (even when
C  "soft failure" has been set "on").
C
C  When a catastrophic failure (IER = 911) is detected in one of
C  the routines in RKSUITE, it calls CHKFL with ASK = .FALSE. and
C  ERROR = .TRUE.  This value of ERROR is SAVEd.
C
C     .. Scalar Arguments ..
      LOGICAL           ASK, ERROR
C     .. Local Scalars ..
      LOGICAL           SAVERR
C     .. Save statement ..
      SAVE              SAVERR
C     .. Data statements ..
      DATA              SAVERR/.FALSE./
C     .. Executable Statements ..
C
      IF (ASK) THEN
         ERROR = SAVERR
      ELSE
         SAVERR = ERROR
      END IF
C
      RETURN
      END
