***************************************************************************
*       10        20        30        40        50        60        70
*23456789012345678901234567890123456789012345678901234567890123456789012
*
       PROGRAM MAIN
       DOUBLE PRECISION Y(100),YP(100),T,TOUT,TEND,STPSZ
       DOUBLE PRECISION RELERR,ABSERR
       INTEGER NEQ,K,IFLAG
*      DOUBLE PRECISION T1,T2,T3(2),T4(2)
       DOUBLE PRECISION KK, CC, MM, FF
       COMMON /PARAMETERS/ KK, CC, MM, FF
       EXTERNAL F
       DATA T,Y,YP / 0.0D0,100*0.0D0,100*0.0D0 /
       DATA TEND,STPSZ,K,NEQ,RELERR,ABSERR / 2*0.0D0,2*0,2*0.0D0 /
       DATA KK, CC, MM, FF / 4*0.0D0 /
*
       OPEN(1,FILE='inputs')
       OPEN(2,FILE='outputs')
       read(1,*) T, K, STPSZ, TOUT, TEND, RELERR, ABSERR
       read(1,*)
       read(1,200) NEQ
  200  format(5x,i10)
       read(1,205) MM
  205  format(5x,1e15.6)
       read(1,210) KK
  210  format(5x,1e15.6)
       read(1,215) CC
  215  format(5x,1e15.6)
       read(1,220) FF
  220  format(5x,1e15.6)
       read(1,225) V0
  225  format(12x,1e15.6)
*      print*, T, K, STPSZ, TOUT, TEND, RELERR, ABSERR
*      print*,'       neq = ', NEQ
*      print*,'        MM = ', MM
*      print*,'        KK = ', KK
*      print*,'        CC = ',CC
*      print*,'        FF = ',FF
*      print*,'speed(t=0) = ',V0
*
* Apply initial conditions
*
       Y(1)  = V0
       Y(2)  = 0.0
       YP(1) = V0
       YP(2) = 0.0
*
* ETIME used for timing cpu execution time
*
*      T1 = ETIME(T3)
*
*      T = 0.0
*      K = 10
*      STPSZ = 0.01
*      TOUT = 0.1
*      TEND = 2.0
*
* IFLAG initializes the integrator
*
       IFLAG = 1
*
* tolerances for integration
*
*      RELERR = 0.00001D0
*      ABSERR = 0.00001D0
*
* write output at time zero
*
       WRITE(2,45) T,Y(1),Y(2)
  45   FORMAT(3E15.6)
*
   10  CONTINUE
*
* control output print frequency
*
       PRINT*,'iflag,tout', IFLAG,T
       TOUT = T + K*STPSZ
*
       CALL DE2(F,NEQ,Y,T,TOUT,RELERR,ABSERR,IFLAG)
*
* write output at each output time interval
*
       WRITE(2,45) T,Y(1),Y(2)
*
*
* If time is left in the simulation take the next time step.
*
       IF (T .GE. TEND) GOTO 100
       GOTO 10
*
* Determine the amount of simulation CPU time using ETIME.
*
* 100  T2 = ETIME(T4)
  100  CONTINUE
       PRINT*,'iflag,tout', IFLAG,T
*      PRINT*,'elapsed cpu time (sec):',T2-T1
*
*
       END
****************************************************************************
       SUBROUTINE F(T,Y,YP)
       DOUBLE PRECISION T,Y(100),YP(100)
       DOUBLE PRECISION PI
       DOUBLE PRECISION X1, X2, X1D
       DOUBLE PRECISION KK, CC, MM, FF
       COMMON /PARAMETERS/ KK, CC, MM, FF
       DATA PI / 3.14159265359D0 /
*
* Update working vectors with integrated state variables
*
       X2 = Y(1)              ! displacements
       X1 = Y(2)              ! velocity
*
* RHS
*
       X1D = -(KK*X2 + CC*X1)/MM + FF/MM
*
* initialize the state vector deriviative
*
       Y(1) = X2             ! translational position
       Y(2) = X1             ! translational velocity
*
       YP(1) = X1             ! translational velocity
       YP(2) = X1D            ! translational acceleration
*
       RETURN
*
       END
*
* de2_iprdct.f 
      SUBROUTINE DE2(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG)               DE     1
C                                                                       DE     2
C   SUBROUTINE  DE  INTEGRATES A SYSTEM OF UP TO 20 FIRST ORDER ORDINARYDE     3
C   DIFFERENTIAL EQUATIONS OF THE FORM                                  DE     4
C             DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN))                     DE     5
C             Y(I) GIVEN AT  T .                                        DE     6
C   THE SUBROUTINE INTEGRATES FROM  T  TO  TOUT .  ON RETURN THE        DE     7
C   PARAMETERS IN THE CALL LIST ARE INITIALIZED FOR CONTINUING THE      DE     8
C   INTEGRATION.  THE USER HAS ONLY TO DEFINE A NEW VALUE  TOUT         DE     9
C   AND CALL  DE  AGAIN.                                                DE    10
C                                                                       DE    11
C   DE  CALLS TWO CODES, THE INTEGRATOR  STEP  AND THE INTERPOLATION    DE    12
C   ROUTINE  INTRP .  STEP  USES A MODIFIED DIVIDED DIFFERENCE FORM OF  DE    13
C   THE ADAMS PECE FORMULAS AND LOCAL EXTRAPOLATION.  IT ADJUSTS THE    DE    14
C   ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR.  NORMALLY EACH CALL DE    15
C   TO  STEP  ADVANCES THE SOLUTION ONE STEP IN THE DIRECTION OF  TOUT .DE    16
C   FOR REASONS OF EFFICIENCY  DE  INTEGRATES BEYOND  TOUT  INTERNALLY, DE    17
C   THOUGH NEVER BEYOND T + 10*(TOUT-T), AND CALLS  INTRP  TO           DE    18
C   INTERPOLATE THE SOLUTION AT  TOUT .  AN OPTION IS PROVIDED TO STOP  DE    19
C   THE INTEGRATION AT  TOUT  BUT IT SHOULD BE USED ONLY IF IT IS       DE    20
C   IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND  TOUT .               DE    21
C                                                                       DE    22
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,       DE    23
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS:  THE INITIAL  DE    24
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.                  DE    25
C                                                                       DE    26
C   THE PARAMETERS FOR  DE  ARE:                                        DE    27
C      F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT DE    28
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     DE    29
C      Y(*) -- SOLUTION VECTOR AT  T                                    DE    30
C      T -- INDEPENDENT VARIABLE                                        DE    31
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED                       DE    32
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCALDE    33
C           ERROR TEST.  AT EACH STEP THE CODE REQUIRES                 DE    34
C             ABS(LOCAL ERROR) .LE. ABS(Y)*RELERR + ABSERR              DE    35
C           FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS  DE    36
C      IFLAG -- INDICATES STATUS OF INTEGRATION                         DE    37
C                                                                       DE    38
C   FIRST CALL TO DE --                                                 DE    39
C                                                                       DE    40
C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAY  DE    41
C   IN THE CALL LIST, Y(NEQN), DECLARE  F  IN AN EXTERNAL STATEMENT,    DE    42
C   SUPPLY THE SUBROUTINE F(T,Y,YP) TO EVALUATE                         DE    43
C      DY(I)/DT = YP(I) = F(T,Y(1),Y(2),...,Y(NEQN))                    DE    44
C   AND INITIALIZE THE PARAMETERS:                                      DE    45
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     DE    46
C      Y(*) -- VECTOR OF INITIAL CONDITIONS                             DE    47
C      T -- STARTING POINT OF INTEGRATION                               DE    48
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED                       DE    49
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES    DE    50
C      IFLAG -- +1,-1.  INDICATOR TO INITIALIZE THE CODE.  NORMAL INPUT DE    51
C           IS +1.  THE USER SHOULD SET IFLAG=-1 ONLY IF IT IS          DE    52
C           IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND  TOUT .       DE    53
C   ALL PARAMETERS EXCEPT  F ,  NEQN  AND  TOUT  MAY BE ALTERED BY THE  DE    54
C   CODE ON OUTPUT SO MUST BE VARIABLES IN THE CALLING PROGRAM.         DE    55
C                                                                       DE    56
C   OUTPUT FROM  DE --                                                  DE    57
C                                                                       DE    58
C      NEQN -- UNCHANGED                                                DE    59
C      Y(*) -- SOLUTION AT  T                                           DE    60
C      T -- LAST POINT REACHED IN INTEGRATION.  NORMAL RETURN HAS       DE    61
C           T = TOUT .                                                  DE    62
C      TOUT -- UNCHANGED                                                DE    63
C      RELERR,ABSERR -- NORMAL RETURN HAS TOLERANCES UNCHANGED.  IFLAG=3DE    64
C           SIGNALS TOLERANCES INCREASED                                DE    65
C      IFLAG = 2 -- NORMAL RETURN.  INTEGRATION REACHED  TOUT           DE    66
C            = 3 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE ERROR      DE    67
C                   TOLERANCES TOO SMALL.  RELERR ,  ABSERR  INCREASED  DE    68
C                   APPROPRIATELY FOR CONTINUING                        DE    69
C            = 4 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE MORE THAN  DE    70
C                   MAXNUM  STEPS NEEDED                                DE    71
C            = 5 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE EQUATIONS  DE    72
C                   APPEAR TO BE STIFF                                  DE    73
C            = 6 -- INVALID INPUT PARAMETERS (FATAL ERROR)              DE    74
C           THE VALUE OF  IFLAG  IS RETURNED NEGATIVE WHEN THE INPUT    DE    75
C           VALUE IS NEGATIVE AND THE INTEGRATION DOES NOT REACH  TOUT ,DE    76
C           I.E., -3, -4, -5.                                           DE    77
C                                                                       DE    78
C   SUBSEQUENT CALLS TO  DE --                                          DE    79
C                                                                       DE    80
C   SUBROUTINE  DE  RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE     DE    81
C   THE INTEGRATION.  IF THE INTEGRATION REACHED  TOUT , THE USER NEED  DE    82
C   ONLY DEFINE A NEW  TOUT  AND CALL AGAIN.  IF THE INTEGRATION DID NOTDE    83
C   REACH  TOUT  AND THE USER WANTS TO CONTINUE, HE JUST CALLS AGAIN.   DE    84
C   THE OUTPUT VALUE OF  IFLAG  IS THE APPROPRIATE INPUT VALUE FOR      DE    85
C   SUBSEQUENT CALLS.  THE ONLY SITUATION IN WHICH IT SHOULD BE ALTERED DE    86
C   IS TO STOP THE INTEGRATION INTERNALLY AT THE NEW  TOUT , I.E.,      DE    87
C   CHANGE OUTPUT  IFLAG=2  TO INPUT  IFLAG=-2 .  ERROR TOLERANCES MAY  DE    88
C   BE CHANGED BY THE USER BEFORE CONTINUING.  ALL OTHER PARAMETERS MUSTDE    89
C   REMAIN UNCHANGED.                                                   DE    90
C                                                                       DE    91
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                        
      LOGICAL START,CRASH,STIFF                                         DE    92
      DIMENSION Y(300),PSI(12)                                          DE    93
      DIMENSION YY(300),WT(300),PHI(300,16),P(300),YP(300),YPOUT(300)   DE    94
      EXTERNAL F                                                        DE    95
      SAVE                                                                      
      COMMON/IPRDCT/IPRDCT
C                                                                       DE    96
C***********************************************************************DE    97
C*  THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT   *DE    98
C*  ROUNDOFF ERROR  U  WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT *DE    99
C*  1.0+U .GT. 1.0 .  U  MUST BE CALCULATED AND  FOURU=4.0*U  INSERTED *DE   100
C*  IN THE FOLLOWING DATA STATEMENT BEFORE USING  DE .  THE ROUTINE    *DE   101
C*  MACHIN  CALCULATES  U .  FOURU  AND  TWOU=2.0*U  MUST ALSO BE      *DE   102
C*  INSERTED IN SUBROUTINE  STEP  BEFORE CALLING  DE .                 *DE   103
      DATA FOURU/8.88D-16/                                              DE   104
C***********************************************************************DE   105
C                                                                       DE   106
C   THE CONSTANT  MAXNUM  IS THE MAXIMUM NUMBER OF STEPS ALLOWED IN ONE DE   107
C   CALL TO  DE .  THE USER MAY CHANGE THIS LIMIT BY ALTERING THE       DE   108
C   FOLLOWING STATEMENT                                                 DE   109
      DATA MAXNUM/500/                                                  DE   110
C   THIS VERSION OF  DE  WILL HANDLE UP TO 20 EQUATIONS.  TO CHANGE THISDE   111
C   NUMBER, ONLY THE NUMBER 20 IN THE DIMENSION STATEMENT ABOVE AND THE DE   112
C   NUMBER 20 IN THE FOLLOWING STATEMENT NEED BE CHANGED                DE   113
C                                                                       DE   114
C            ***            ***            ***                          DE   115
C   TEST FOR IMPROPER PARAMETERS                                        DE   116
C                                                                       DE   117
      IF(NEQN .LT. 1  .OR.  NEQN .GT. 300) GO TO 10                     DE   118
      IF(T .EQ. TOUT) GO TO 10                                          DE   119
      IF(RELERR .LT. 0.D0 .OR.  ABSERR .LT. 0.D0)GO TO 10               DE   120
      EPS = DMAX1(RELERR,ABSERR)                                        DE   121
      IF(EPS .LE. 0.D0)GO TO 10                                         DE   122
      IF(IFLAG .EQ. 0) GO TO 10                                         DE   123
      ISN = ISIGN(1,IFLAG)                                              DE   124
      IFLAG = IABS(IFLAG)                                               DE   125
      IF(IFLAG .EQ. 1) GO TO 20                                         DE   126
      IF(T .NE. TOLD) GO TO 10                                          DE   127
      IF(IFLAG .GE. 2  .AND.  IFLAG .LE. 5) GO TO 20                    DE   128
   10 IFLAG = 6                                                         DE   129
      RETURN                                                            DE   130
C                                                                       DE   131
C   ON EACH CALL SET INTERVAL OF INTEGRATION AND COUNTER FOR NUMBER OF  DE   132
C   STEPS.  ADJUST INPUT ERROR TOLERANCES TO DEFINE WEIGHT VECTOR FOR   DE   133
C   SUBROUTINE  STEP                                                    DE   134
C                                                                       DE   135
   20 DEL = TOUT - T                                                    DE   136
      ABSDEL =DABS(DEL)                                                 DE   137
      TEND = T + 10.0D0*DEL                                             DE   138
      IF(ISN .LT. 0) TEND = TOUT                                        DE   139
      NOSTEP = 0                                                        DE   140
      KLE4 = 0                                                          DE   141
      STIFF = .FALSE.                                                   DE   142
      RELEPS = RELERR/EPS                                               DE   143
      ABSEPS = ABSERR/EPS                                               DE   144
      IF(IFLAG .EQ. 1) GO TO 30                                         DE   145
      IF(ISNOLD .LT. 0) GO TO 30                                        DE   146
      IF(DELSGN*DEL .GT. 0.D0)GO TO 50                                  DE   147
C                                                                       DE   148
C   ON START AND RESTART ALSO SET WORK VARIABLES X AND YY(*), STORE THE DE   149
C   DIRECTION OF INTEGRATION AND INITIALIZE THE STEP SIZE               DE   150
C                                                                       DE   151
   30 START = .TRUE.                                                    DE   152
      X = T                                                             DE   153
      DO 40 L = 1,NEQN                                                  DE   154
   40   YY(L) = Y(L)                                                    DE   155
      DELSGN =DSIGN(1.D0,DEL)                                           DE   156
      H=DSIGN(DMAX1(DABS(TOUT-X),FOURU*DABS(X)),TOUT-X)                 DE   157
C                                                                       DE   158
C   IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN                DE   159
C                                                                       DE   160
   50 IF(DABS(X-T).LT. ABSDEL) GO TO 60                                 DE   161
      CALL INTRP2(X,YY,TOUT,Y,YPOUT,NEQN,KOLD,PHI,PSI)                   DE   162
      IFLAG = 2                                                         DE   163
      T = TOUT                                                          DE   164
      TOLD = T                                                          DE   165
      ISNOLD = ISN                                                      DE   166
      RETURN                                                            DE   167
C                                                                       DE   168
C   IF CANNOT GO PAST OUTPUT POINT AND SUFFICIENTLY CLOSE,              DE   169
C   EXTRAPOLATE AND RETURN                                              DE   170
C                                                                       DE   171
   60 IF(ISN .GT. 0  .OR. DABS(TOUT-X) .GE. FOURU*DABS(X))GO TO 80      DE   172
      H = TOUT - X                                                      DE   173
      IPRDCT=0
      CALL F(X,YY,YP)                                                   DE   174
      DO 70 L = 1,NEQN                                                  DE   175
   70   Y(L) = YY(L) + H*YP(L)                                          DE   176
      IFLAG = 2                                                         DE   177
      T = TOUT                                                          DE   178
      TOLD = T                                                          DE   179
      ISNOLD = ISN                                                      DE   180
      RETURN                                                            DE   181
C                                                                       DE   182
C   TEST FOR TOO MUCH WORK                                              DE   183
C                                                                       DE   184
   80 IF(NOSTEP .LT. MAXNUM) GO TO 100                                  DE   185
      IFLAG = ISN*4                                                     DE   186
      IF(STIFF) IFLAG = ISN*5                                           DE   187
      DO 90 L = 1,NEQN                                                  DE   188
   90   Y(L) = YY(L)                                                    DE   189
      T = X                                                             DE   190
      TOLD = T                                                          DE   191
      ISNOLD = 1                                                        DE   192
      RETURN                                                            DE   193
C                                                                       DE   194
C   LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP                  DE   195
C                                                                       DE   196
  100 H =DSIGN(DMIN1(DABS(H),DABS(TEND-X)),H)                           DE   197
      DO 110 L = 1,NEQN                                                 DE   198
  110   WT(L) =RELEPS*DABS(YY(L)) + ABSEPS                              DE   199
      CALL STEPOCP(X,YY,F,NEQN,H,EPS,WT,START,                             DE   200
     1  HOLD,K,KOLD,CRASH,PHI,P,YP,PSI)                                 DE   201
C                                                                       DE   202
C   TEST FOR TOLERANCES TOO SMALL                                       DE   203
C                                                                       DE   204
      IF(.NOT.CRASH) GO TO 130                                          DE   205
      IFLAG = ISN*3                                                     DE   206
      RELERR = EPS*RELEPS                                               DE   207
      ABSERR = EPS*ABSEPS                                               DE   208
      DO 120 L = 1,NEQN                                                 DE   209
  120   Y(L) = YY(L)                                                    DE   210
      T = X                                                             DE   211
      TOLD = T                                                          DE   212
      ISNOLD = 1                                                        DE   213
      RETURN                                                            DE   214
C                                                                       DE   215
C   AUGMENT COUNTER ON WORK AND TEST FOR STIFFNESS                      DE   216
C                                                                       DE   217
  130 NOSTEP = NOSTEP + 1                                               DE   218
      KLE4 = KLE4 + 1                                                   DE   219
      IF(KOLD .GT. 4) KLE4 = 0                                          DE   220
      IF(KLE4 .GE. 50) STIFF = .TRUE.                                   DE   221
      GO TO 50                                                          DE   222
      END                                                               DE   223
      SUBROUTINE STEPOCP(X,Y,F,NEQN,H,EPS,WT,START,                        STEP   1
     1  HOLD,K,KOLD,CRASH,PHI,P,YP,PSI)                                 STEP   2
C                                                                       STEP   3
C   SUBROUTINE  STEP  INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY       STEP   4
C   DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A    STEP   5
C   MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS.  LOCAL STEP   6
C   EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY.   STEP   7
C   THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR STEP   8
c   per unit step in a generalized sense.  special devices are included step   9
C   TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING STEP  10
C   TOO MUCH ACCURACY.                                                  STEP  11
C                                                                       STEP  12
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,       STEP  13
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS:  THE INITIAL  STEP  14
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.                  STEP  15
C                                                                       STEP  16
C                                                                       STEP  17
C   THE PARAMETERS REPRESENT:                                           STEP  18
C      X -- INDEPENDENT VARIABLE                                        STEP  19
C      Y(*) -- SOLUTION VECTOR AT X                                     STEP  20
C      YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT  X  AFTER SUCCESSFUL   STEP  21
C           STEP                                                        STEP  22
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     STEP  23
C      H -- APPROPRIATE STEP SIZE FOR NEXT STEP.  NORMALLY DETERMINED BYSTEP  24
C           CODE                                                        STEP  25
C      EPS -- LOCAL ERROR TOLERANCE.  MUST BE VARIABLE                  STEP  26
C      WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION                   STEP  27
C      START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP,  .FALSE.    STEP  28
C           OTHERWISE                                                   STEP  29
C      HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP                  STEP  30
C      K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE)        STEP  31
C      KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP                      STEP  32
C      CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN,  STEP  33
C           .FALSE. OTHERWISE.                                          STEP  34
C   THE ARRAYS  PHI, PSI  ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE STEP  35
C   INTRP .  THE ARRAY  P  IS INTERNAL TO THE CODE.                     STEP  36
C                                                                       STEP  37
C   INPUT TO  STEP                                                      STEP  38
C                                                                       STEP  39
C      FIRST CALL --                                                    STEP  40
C                                                                       STEP  41
C   THE USER MUST PROVIDE STORAGE IN HIS DRIVER PROGRAM FOR ALL ARRAYS  STEP  42
C   IN THE CALL LIST, NAMELY                                            STEP  43
C                                                                       STEP  44
C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12)  STEP  45
C                                                                       STEP  46
C   THE USER MUST ALSO DECLARE  START  AND  CRASH  LOGICAL VARIABLES    STEP  47
C   AND  F  AN EXTERNAL SUBROUTINE, SUPPLY THE SUBROUTINE  F(X,Y,YP)    STEP  48
C   TO EVALUATE                                                         STEP  49
C      DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN))                    STEP  50
C   AND INITIALIZE ONLY THE FOLLOWING PARAMETERS:                       STEP  51
C      X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE                   STEP  52
C      Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES          STEP  53
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     STEP  54
C      H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION       STEP  55
C           AND MAXIMUM SIZE OF STEP.  MUST BE VARIABLE                 STEP  56
C      EPS -- LOCAL ERROR TOLERANCE PER STEP.  MUST BE VARIABLE         STEP  57
C      WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION          STEP  58
C      START -- .TRUE.                                                  STEP  59
C                                                                       STEP  60
C   STEP  REQUIRES THE L2 NORM OF THE VECTOR WITH COMPONENTS            STEP  61
C   LOCAL ERROR(L)/WT(L)  BE LESS THAN  EPS  FOR A SUCCESSFUL STEP.  THESTEP  62
C   ARRAY  WT  ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE     STEP  63
C   FOR HIS PROBLEM.  FOR EXAMPLE,                                      STEP  64
C      WT(L) = 1.0  SPECIFIES ABSOLUTE ERROR,                           STEP  65
C            = ABS(Y(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF THESTEP  66
C                 L-TH COMPONENT OF THE SOLUTION,                       STEP  67
C            = ABS(YP(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF   STEP  68
C                 THE L-TH COMPONENT OF THE DERIVATIVE,                 STEP  69
C            = AMAX1(WT(L),ABS(Y(L)))  ERROR RELATIVE TO THE LARGEST    STEP  70
C                 MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR,          STEP  71
C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  SPECIFIES A MIXED     STEP  72
C                 RELATIVE-ABSOLUTE TEST WHERE  RELERR  IS RELATIVE     STEP  73
C                 ERROR,  ABSERR  IS ABSOLUTE ERROR AND  EPS =          STEP  74
C                 AMAX1(RELERR,ABSERR) .                                STEP  75
C                                                                       STEP  76
C      SUBSEQUENT CALLS --                                              STEP  77
C                                                                       STEP  78
C   SUBROUTINE  STEP  IS DESIGNED SO THAT ALL INFORMATION NEEDED TO     STEP  79
C   CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE  H  AND THE ORDER STEP  80
C   K , IS RETURNED WITH EACH STEP.  WITH THE EXCEPTION OF THE STEP     STEP  81
C   SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS  STEP  82
C   SHOULD BE ALTERED.  THE ARRAY  WT  MUST BE UPDATED AFTER EACH STEP  STEP  83
C   TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE.  NORMALLY THE    STEP  84
C   INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE   STEP  85
C   SOLUTION INTERPOLATED THERE WITH SUBROUTINE  INTRP .  IF IT IS      STEP  86
C   IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE   STEP  87
C   REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP     STEP  88
C   LARGER THAN THE  H  INPUT.  CHANGING THE DIRECTION OF INTEGRATION,  STEP  89
C   I.E., THE SIGN OF  H , REQUIRES THE USER SET  START = .TRUE. BEFORE STEP  90
C   CALLING  STEP  AGAIN.  THIS IS THE ONLY SITUATION IN WHICH  START   STEP  91
C   SHOULD BE ALTERED.                                                  STEP  92
C                                                                       STEP  93
C   OUTPUT FROM  STEP                                                   STEP  94
C                                                                       STEP  95
C      SUCCESSFUL STEP --                                               STEP  96
C                                                                       STEP  97
C   THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH  START  AND  STEP  98
C   CRASH  SET .FALSE. .  X  REPRESENTS THE INDEPENDENT VARIABLE        STEP  99
C   ADVANCED ONE STEP OF LENGTH  HOLD  FROM ITS VALUE ON INPUT AND  Y   STEP 100
C   THE SOLUTION VECTOR AT THE NEW VALUE OF  X .  ALL OTHER PARAMETERS  STEP 101
C   REPRESENT INFORMATION CORRESPONDING TO THE NEW  X  NEEDED TO        STEP 102
C   CONTINUE THE INTEGRATION.                                           STEP 103
C                                                                       STEP 104
C      UNSUCCESSFUL STEP --                                             STEP 105
C                                                                       STEP 106
C   WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION,    STEP 107
C   THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND  CRASH = .TRUE. .  STEP 108
C   AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE     STEP 109
C   ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT       STEP 110
C   BEFORE RETURNING.  TO CONTINUE WITH THE LARGER TOLERANCE, THE USER  STEP 111
C   JUST CALLS THE CODE AGAIN.  A RESTART IS NEITHER REQUIRED NOR       STEP 112
C   DESIRABLE.                                                          STEP 113
C                                                                       STEP 114
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                        
      LOGICAL START,CRASH,PHASE1,NORND                                  STEP 115
      DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12)  STEP 116
      DIMENSION ALPHA(12),BETA(12),SIG(13),W(12),V(12),G(13),           STEP 117
     1  GSTR(13),TWO(13)                                                STEP 118
      COMMON/IPRDCT/IPRDCT
      EXTERNAL F                                                        STEP 119
      SAVE                                                                      
C***********************************************************************STEP 120
C*  THE ONLY MACHINE DEPENDENT CONSTANTS ARE BASED ON THE MACHINE UNIT *STEP 121
C*  ROUNDOFF ERROR  U  WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT *STEP 122
C*  1.0+U .GT. 1.0  .  THE USER MUST CALCULATE  U  AND INSERT          *STEP 123
C*  TWOU=2.0*U  AND  FOURU=4.0*U  IN THE DATA STATEMENT BEFORE CALLING *STEP 124
C*  THE CODE.  THE ROUTINE  MACHIN  CALCULATES  U .                    *STEP 125
      DATA TWOU,FOURU/2.8D-14,5.68D-14/                                 STEP 126
C***********************************************************************STEP 127
      DATA TWO/2D0,4D0,8D0,16D0,32D0,64D0,128D0,256D0,512D0,1024D0,     STEP 128
     1  2048D0,4096D0,8192D0/                                           STEP 129
      DATA GSTR/.5D0,.0833D0,.0417D0,.0264D0,.0188D0,.0143D0,.0114D0,   STEP 130
     1 .936D-2,.789D-2,.679D-2,.592D-2,.524D-2,.468D-2/                 STEP 131
      DATA G(1),G(2)/1.0D0,0.5D0/,SIG(1)/1.0D0/                         STEP 132
C                                                                       STEP 133
C                                                                       STEP 134
C       ***     BEGIN BLOCK 0     ***                                   STEP 135
C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE      STEP 136
C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A      STEP 137
C   STARTING STEP SIZE.                                                 STEP 138
C                   ***                                                 STEP 139
C                                                                       STEP 140
C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE              STEP 141
C                                                                       STEP 142
      CRASH = .TRUE.                                                    STEP 143
      IF(DABS(H) .GE. FOURU*DABS(X)) GO TO 5                            STEP 144
      H =DSIGN(FOURU*DABS(X),H)                                         STEP 145
      RETURN                                                            STEP 146
    5 P5EPS = 0.5D0*EPS                                                 STEP 147
C                                                                       STEP 148
C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE STEP 149
C                                                                       STEP 150
      ROUND = 0.0D0                                                     STEP 151
      DO 10 L = 1,NEQN                                                  STEP 152
   10   ROUND = ROUND + (Y(L)/WT(L))**2                                 STEP 153
      ROUND =TWOU*DSQRT(ROUND)                                          STEP 154
      IF(P5EPS .GE. ROUND) GO TO 15                                     STEP 155
      EPS = 2D0*ROUND*(1D0 + FOURU)                                     STEP 156
      RETURN                                                            STEP 157
   15 CRASH = .FALSE.                                                   STEP 158
      IF(.NOT.START) GO TO 99                                           STEP 159
C                                                                       STEP 160
C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP           STEP 161
C                                                                       STEP 162
      IPRDCT=0
      CALL F(X,Y,YP)                                                    STEP 163
      SUM = 0.0D0                                                       STEP 164
      DO 20 L = 1,NEQN                                                  STEP 165
        PHI(L,1) = YP(L)                                                STEP 166
        PHI(L,2) = 0.0D0                                                STEP 167
   20   SUM = SUM + (YP(L)/WT(L))**2                                    STEP 168
      SUM =DSQRT(SUM)                                                   STEP 169
      ABSH =DABS(H)                                                     STEP 170
      IF(EPS .LT. 16D0*SUM*H*H) ABSH = 0.25D0*DSQRT(EPS/SUM)            STEP 171
      H =DSIGN(DMAX1(ABSH,FOURU*DABS(X)),H)                             STEP 172
      HOLD = 0.0D0                                                      STEP 173
      K = 1                                                             STEP 174
      KOLD = 0                                                          STEP 175
      START = .FALSE.                                                   STEP 176
      PHASE1 = .TRUE.                                                   STEP 177
      NORND = .TRUE.                                                    STEP 178
      IF(P5EPS .GT. 100D0*ROUND) GO TO 99                               STEP 179
      NORND = .FALSE.                                                   STEP 180
      DO 25 L = 1,NEQN                                                  STEP 181
   25   PHI(L,15) = 0.0D0                                               STEP 182
   99 IFAIL = 0                                                         STEP 183
C       ***     END BLOCK 0     ***                                     STEP 184
C                                                                       STEP 185
C       ***     BEGIN BLOCK 1     ***                                   STEP 186
C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING    STEP 187
C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.         STEP 188
C                   ***                                                 STEP 189
C                                                                       STEP 190
  100 KP1 = K+1                                                         STEP 191
      KP2 = K+2                                                         STEP 192
      KM1 = K-1                                                         STEP 193
      KM2 = K-2                                                         STEP 194
C                                                                       STEP 195
C   NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT  STEP 196
C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE                          STEP 197
C                                                                       STEP 198
      IF(H .NE. HOLD) NS = 0                                            STEP 199
      IF(NS .LE. KOLD) NS=NS+1                                          STEP 200
      NSP1 = NS+1                                                       STEP 201
      IF (K .LT. NS) GO TO 199                                          STEP 202
C                                                                       STEP 203
C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH    STEP 204
C   ARE CHANGED                                                         STEP 205
C                                                                       STEP 206
      BETA(NS) = 1.0D0                                                  STEP 207
      REALNS = NS                                                       STEP 208
      ALPHA(NS) =1.D0/REALNS                                            STEP 209
      TEMP1 = H*REALNS                                                  STEP 210
      SIG(NSP1) = 1.0D0                                                 STEP 211
      IF(K .LT. NSP1) GO TO 110                                         STEP 212
      DO 105 I = NSP1,K                                                 STEP 213
        IM1 = I-1                                                       STEP 214
        TEMP2 = PSI(IM1)                                                STEP 215
        PSI(IM1) = TEMP1                                                STEP 216
        BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2                              STEP 217
        TEMP1 = TEMP2 + H                                               STEP 218
        ALPHA(I) = H/TEMP1                                              STEP 219
        REALI = I                                                       STEP 220
  105   SIG(I+1) = REALI*ALPHA(I)*SIG(I)                                STEP 221
  110 PSI(K) = TEMP1                                                    STEP 222
C                                                                       STEP 223
C   COMPUTE COEFFICIENTS G(*)                                           STEP 224
C                                                                       STEP 225
C   INITIALIZE V(*) AND SET W(*).  G(2) IS SET IN DATA STATEMENT        STEP 226
C                                                                       STEP 227
      IF(NS .GT. 1) GO TO 120                                           STEP 228
      DO 115 IQ = 1,K                                                   STEP 229
        TEMP3 = IQ*(IQ+1)                                               STEP 230
        V(IQ) =1.D0/TEMP3                                               STEP 231
  115   W(IQ) = V(IQ)                                                   STEP 232
      GO TO 140                                                         STEP 233
C                                                                       STEP 234
C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)                   STEP 235
C                                                                       STEP 236
  120 IF(K .LE. KOLD) GO TO 130                                         STEP 237
      TEMP4 = K*KP1                                                     STEP 238
      V(K) =1.D0/TEMP4                                                  STEP 239
      NSM2 = NS-2                                                       STEP 240
      IF(NSM2 .LT. 1) GO TO 130                                         STEP 241
      DO 125 J = 1,NSM2                                                 STEP 242
        I = K-J                                                         STEP 243
  125   V(I) = V(I) - ALPHA(J+1)*V(I+1)                                 STEP 244
C                                                                       STEP 245
C   UPDATE V(*) AND SET W(*)                                            STEP 246
C                                                                       STEP 247
  130 LIMIT1 = KP1 - NS                                                 STEP 248
      TEMP5 = ALPHA(NS)                                                 STEP 249
      DO 135 IQ = 1,LIMIT1                                              STEP 250
        V(IQ) = V(IQ) - TEMP5*V(IQ+1)                                   STEP 251
  135   W(IQ) = V(IQ)                                                   STEP 252
      G(NSP1) = W(1)                                                    STEP 253
C                                                                       STEP 254
C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)                            STEP 255
C                                                                       STEP 256
  140 NSP2 = NS + 2                                                     STEP 257
      IF(KP1 .LT. NSP2) GO TO 199                                       STEP 258
      DO 150 I = NSP2,KP1                                               STEP 259
        LIMIT2 = KP2 - I                                                STEP 260
        TEMP6 = ALPHA(I-1)                                              STEP 261
        DO 145 IQ = 1,LIMIT2                                            STEP 262
  145     W(IQ) = W(IQ) - TEMP6*W(IQ+1)                                 STEP 263
  150   G(I) = W(1)                                                     STEP 264
  199   CONTINUE                                                        STEP 265
C       ***     END BLOCK 1     ***                                     STEP 266
C                                                                       STEP 267
C       ***     BEGIN BLOCK 2     ***                                   STEP 268
C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED       STEP 269
C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,   STEP 270
C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.                        STEP 271
C                   ***                                                 STEP 272
C                                                                       STEP 273
C   CHANGE PHI TO PHI STAR                                              STEP 274
C                                                                       STEP 275
      IF(K .LT. NSP1) GO TO 215                                         STEP 276
      DO 210 I = NSP1,K                                                 STEP 277
        TEMP1 = BETA(I)                                                 STEP 278
        DO 205 L = 1,NEQN                                               STEP 279
  205     PHI(L,I) = TEMP1*PHI(L,I)                                     STEP 280
  210   CONTINUE                                                        STEP 281
C                                                                       STEP 282
C   PREDICT SOLUTION AND DIFFERENCES                                    STEP 283
C                                                                       STEP 284
  215 DO 220 L = 1,NEQN                                                 STEP 285
        PHI(L,KP2) = PHI(L,KP1)                                         STEP 286
        PHI(L,KP1) = 0.0D0                                              STEP 287
  220   P(L) = 0.0D0                                                    STEP 288
      DO 230 J = 1,K                                                    STEP 289
        I = KP1 - J                                                     STEP 290
        IP1 = I+1                                                       STEP 291
        TEMP2 = G(I)                                                    STEP 292
        DO 225 L = 1,NEQN                                               STEP 293
          P(L) = P(L) + TEMP2*PHI(L,I)                                  STEP 294
  225     PHI(L,I) = PHI(L,I) + PHI(L,IP1)                              STEP 295
  230   CONTINUE                                                        STEP 296
      IF(NORND) GO TO 240                                               STEP 297
      DO 235 L = 1,NEQN                                                 STEP 298
        TAU = H*P(L) - PHI(L,15)                                        STEP 299
        P(L) = Y(L) + TAU                                               STEP 300
  235   PHI(L,16) = (P(L) - Y(L)) - TAU                                 STEP 301
      GO TO 250                                                         STEP 302
  240 DO 245 L = 1,NEQN                                                 STEP 303
  245   P(L) = Y(L) + H*P(L)                                            STEP 304
  250 XOLD = X                                                          STEP 305
      X = X + H                                                         STEP 306
      ABSH =DABS(H)                                                     STEP 307
      IPRDCT=IPRDCT+1
      CALL F(X,P,YP)                                                    STEP 308
C                                                                       STEP 309
C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2                                 STEP 310
C                                                                       STEP 311
      ERKM2 = 0.0D0                                                     STEP 312
      ERKM1 = 0.0D0                                                     STEP 313
      ERK = 0.0D0                                                       STEP 314
      DO 265 L = 1,NEQN                                                 STEP 315
        TEMP3 =1.D0/WT(L)                                               STEP 316
        TEMP4 = YP(L) - PHI(L,1)                                        STEP 317
        IF(KM2)265,260,255                                              STEP 318
  255   ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2                   STEP 319
  260   ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2                     STEP 320
  265   ERK = ERK + (TEMP4*TEMP3)**2                                    STEP 321
      IF(KM2)280,275,270                                                STEP 322
  270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*DSQRT(ERKM2)                      STEP 323
  275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*DSQRT(ERKM1)                        STEP 324
  280 TEMP5 = ABSH*DSQRT(ERK)                                           STEP 325
      ERR = TEMP5*(G(K)-G(KP1))                                         STEP 326
      ERK = TEMP5*SIG(KP1)*GSTR(K)                                      STEP 327
      KNEW = K                                                          STEP 328
C                                                                       STEP 329
C   TEST IF ORDER SHOULD BE LOWERED                                     STEP 330
C                                                                       STEP 331
      IF(KM2)299,290,285                                                STEP 332
  285 IF(DMAX1(ERKM1,ERKM2) .LE. ERK) KNEW = KM1                        STEP 333
      GO TO 299                                                         STEP 334
  290 IF(ERKM1 .LE. .5D0*ERK)KNEW = KM1                                 STEP 335
C                                                                       STEP 336
C   TEST IF STEP SUCCESSFUL                                             STEP 337
C                                                                       STEP 338
  299 IF(ERR .LE. EPS) GO TO 400                                        STEP 339
C       ***     END BLOCK 2     ***                                     STEP 340
C                                                                       STEP 341
C       ***     BEGIN BLOCK 3     ***                                   STEP 342
C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .           STEP 343
C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE STEP 344
C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR      STEP 345
C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINESTEP 346
C   PRECISION.                                                          STEP 347
C                   ***                                                 STEP 348
C                                                                       STEP 349
C   RESTORE X, PHI(*,*) AND PSI(*)                                      STEP 350
C                                                                       STEP 351
      PHASE1 = .FALSE.                                                  STEP 352
      X = XOLD                                                          STEP 353
      DO 310 I = 1,K                                                    STEP 354
        TEMP1 =1.D0/BETA(I)                                             STEP 355
        IP1 = I+1                                                       STEP 356
        DO 305 L = 1,NEQN                                               STEP 357
  305     PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))                      STEP 358
  310   CONTINUE                                                        STEP 359
      IF(K .LT. 2) GO TO 320                                            STEP 360
      DO 315 I = 2,K                                                    STEP 361
  315   PSI(I-1) = PSI(I) - H                                           STEP 362
C                                                                       STEP 363
C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP   STEP 364
C   SIZE                                                                STEP 365
C                                                                       STEP 366
  320 IFAIL = IFAIL + 1                                                 STEP 367
      TEMP2 = 0.5D0                                                     STEP 368
      IF(IFAIL - 3) 335,330,325                                         STEP 369
  325 IF(P5EPS .LT. .25D0*ERK)TEMP2 =DSQRT(P5EPS/ERK)                   STEP 370
  330 KNEW = 1                                                          STEP 371
  335 H = TEMP2*H                                                       STEP 372
      K = KNEW                                                          STEP 373
      IF(DABS(H) .GE. FOURU*DABS(X)) GO TO 340                          STEP 374
      CRASH = .TRUE.                                                    STEP 375
      H =DSIGN(FOURU*DABS(X),H)                                         STEP 376
      EPS = EPS + EPS                                                   STEP 377
      RETURN                                                            STEP 378
  340 GO TO 100                                                         STEP 379
C       ***     END BLOCK 3     ***                                     STEP 380
C                                                                       STEP 381
C       ***     BEGIN BLOCK 4     ***                                   STEP 382
C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE   STEP 383
C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE         STEP 384
C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.     STEP 385
C                   ***                                                 STEP 386
  400 KOLD = K                                                          STEP 387
      HOLD = H                                                          STEP 388
C                                                                       STEP 389
C   CORRECT AND EVALUATE                                                STEP 390
C                                                                       STEP 391
      TEMP1 = H*G(KP1)                                                  STEP 392
      IF(NORND) GO TO 410                                               STEP 393
      DO 405 L = 1,NEQN                                                 STEP 394
        RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)                      STEP 395
        Y(L) = P(L) + RHO                                               STEP 396
  405   PHI(L,15) = (Y(L) - P(L)) - RHO                                 STEP 397
      GO TO 420                                                         STEP 398
  410 DO 415 L = 1,NEQN                                                 STEP 399
  415   Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))                          STEP 400
  420 CONTINUE
      IPRDCT=0
      CALL F(X,Y,YP)                                                    STEP 401
C                                                                       STEP 402
C   UPDATE DIFFERENCES FOR NEXT STEP                                    STEP 403
C                                                                       STEP 404
      DO 425 L = 1,NEQN                                                 STEP 405
        PHI(L,KP1) = YP(L) - PHI(L,1)                                   STEP 406
  425   PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)                            STEP 407
      DO 435 I = 1,K                                                    STEP 408
        DO 430 L = 1,NEQN                                               STEP 409
  430     PHI(L,I) = PHI(L,I) + PHI(L,KP1)                              STEP 410
  435   CONTINUE                                                        STEP 411
C                                                                       STEP 412
C   ESTIMATE ERROR AT ORDER K+1 UNLESS:                                 STEP 413
C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,                           STEP 414
C     ALREADY DECIDED TO LOWER ORDER,                                   STEP 415
C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE                     STEP 416
C                                                                       STEP 417
      ERKP1 = 0.0D0                                                     STEP 418
      IF(KNEW .EQ. KM1  .OR.  K .EQ. 12) PHASE1 = .FALSE.               STEP 419
      IF(PHASE1) GO TO 450                                              STEP 420
      IF(KNEW .EQ. KM1) GO TO 455                                       STEP 421
      IF(KP1 .GT. NS) GO TO 460                                         STEP 422
      DO 440 L = 1,NEQN                                                 STEP 423
  440   ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2                           STEP 424
      ERKP1 = ABSH*GSTR(KP1)*DSQRT(ERKP1)                               STEP 425
C                                                                       STEP 426
C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER     STEP 427
C   FOR NEXT STEP                                                       STEP 428
C                                                                       STEP 429
      IF(K .GT. 1) GO TO 445                                            STEP 430
      IF(ERKP1 .GE. .5D0*ERK)GO TO 460                                  STEP 431
      GO TO 450                                                         STEP 432
  445 IF(ERKM1 .LE. DMIN1(ERK,ERKP1)) GO TO 455                         STEP 433
      IF(ERKP1 .GE. ERK  .OR.  K .EQ. 12) GO TO 460                     STEP 434
C                                                                       STEP 435
C   HERE ERKP1 .LT. ERK .LT. AMAX1(ERKM1,ERKM2) ELSE ORDER WOULD HAVE   STEP 436
C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED                STEP 437
C                                                                       STEP 438
C   RAISE ORDER                                                         STEP 439
C                                                                       STEP 440
  450 K = KP1                                                           STEP 441
      ERK = ERKP1                                                       STEP 442
      GO TO 460                                                         STEP 443
C                                                                       STEP 444
C   LOWER ORDER                                                         STEP 445
C                                                                       STEP 446
  455 K = KM1                                                           STEP 447
      ERK = ERKM1                                                       STEP 448
C                                                                       STEP 449
C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP        STEP 450
C                                                                       STEP 451
  460 HNEW = H + H                                                      STEP 452
      IF(PHASE1) GO TO 465                                              STEP 453
      IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465                             STEP 454
      HNEW = H                                                          STEP 455
      IF(P5EPS .GE. ERK) GO TO 465                                      STEP 456
      TEMP2 = K+1                                                       STEP 457
      R = (P5EPS/ERK)**(1D0/TEMP2)                                      STEP 458
      HNEW=ABSH*DMAX1(0.5D0,DMIN1(0.9D0,R))                             STEP 459
      HNEW =DSIGN(DMAX1(HNEW,FOURU*DABS(X)),H)                          STEP 460
  465 H = HNEW                                                          STEP 461
      RETURN                                                            STEP 462
C       ***     END BLOCK 4     ***                                     STEP 463
      END                                                               STEP 464
      SUBROUTINE INTRP2(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI)           INTR   1
C                                                                       INTR   2
C   THE METHODS IN SUBROUTINE  STEP  APPROXIMATE THE SOLUTION NEAR  X   INTR   3
C   BY A POLYNOMIAL.  SUBROUTINE  INTRP  APPROXIMATES THE SOLUTION AT   INTR   4
C   XOUT  BY EVALUATING THE POLYNOMIAL THERE.  INFORMATION DEFINING THISINTR   5
C   POLYNOMIAL IS PASSED FROM  STEP  SO  INTRP  CANNOT BE USED ALONE.   INTR   6
C                                                                       INTR   7
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,       INTR   8
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS:  THE INITIAL  INTR   9
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.                  INTR  10
C                                                                       INTR  11
C   INPUT TO INTRP --                                                   INTR  12
C                                                                       INTR  13
C   THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN  INTR  14
C   THE CALL LIST                                                       INTR  15
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                        
       DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12)    INTR  16
C   AND DEFINES                                                         INTR  17
C      XOUT -- POINT AT WHICH SOLUTION IS DESIRED.                      INTR  18
C   THE REMAINING PARAMETERS ARE DEFINED IN  STEP  AND PASSED TO  INTRP INTR  19
C   FROM THAT SUBROUTINE                                                INTR  20
C                                                                       INTR  21
C   OUTPUT FROM  INTRP --                                               INTR  22
C                                                                       INTR  23
C      YOUT(*) -- SOLUTION AT  XOUT                                     INTR  24
C      YPOUT(*) -- DERIVATIVE OF SOLUTION AT  XOUT                      INTR  25
C   THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT    INTR  26
C   VALUES.  INTEGRATION WITH  STEP  MAY BE CONTINUED.                  INTR  27
C                                                                       INTR  28
      DIMENSION G(13),W(13),RHO(13)                                     INTR  29
      SAVE                                                                      
      DATA G(1)/1D0/,RHO(1)/1D0/                                        INTR  30
C                                                                       INTR  31
      HI = XOUT - X                                                     INTR  32
      KI = KOLD + 1                                                     INTR  33
      KIP1 = KI + 1                                                     INTR  34
C                                                                       INTR  35
C   INITIALIZE W(*) FOR COMPUTING G(*)                                  INTR  36
C                                                                       INTR  37
      DO 5 I = 1,KI                                                     INTR  38
        TEMP1 = I                                                       INTR  39
    5   W(I) =1.D0/TEMP1                                                INTR  40
      TERM = 0.0D0                                                      INTR  41
C                                                                       INTR  42
C   COMPUTE G(*)                                                        INTR  43
C                                                                       INTR  44
      DO 15 J = 2,KI                                                    INTR  45
        JM1 = J - 1                                                     INTR  46
        PSIJM1 = PSI(JM1)                                               INTR  47
        GAMMA = (HI + TERM)/PSIJM1                                      INTR  48
        ETA = HI/PSIJM1                                                 INTR  49
        LIMIT1 = KIP1 - J                                               INTR  50
        DO 10 I = 1,LIMIT1                                              INTR  51
   10     W(I) = GAMMA*W(I) - ETA*W(I+1)                                INTR  52
        G(J) = W(1)                                                     INTR  53
        RHO(J) = GAMMA*RHO(JM1)                                         INTR  54
   15   TERM = PSIJM1                                                   INTR  55
C                                                                       INTR  56
C   INTERPOLATE                                                         INTR  57
C                                                                       INTR  58
      DO 20 L = 1,NEQN                                                  INTR  59
        YPOUT(L) = 0.0D0                                                INTR  60
   20   YOUT(L) = 0.0D0                                                 INTR  61
      DO 30 J = 1,KI                                                    INTR  62
        I = KIP1 - J                                                    INTR  63
        TEMP2 = G(I)                                                    INTR  64
        TEMP3 = RHO(I)                                                  INTR  65
        DO 25 L = 1,NEQN                                                INTR  66
          YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)                            INTR  67
   25     YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)                          INTR  68
   30   CONTINUE                                                        INTR  69
      DO 35 L = 1,NEQN                                                  INTR  70
   35   YOUT(L) = Y(L) + HI*YOUT(L)                                     INTR  71
      RETURN                                                            INTR  72
      END                                                               INTR  73
