      SUBROUTINE RKF45(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK,  RKF 0010
     1           R,D)                                                   RKF 0020
C     FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD                    RKF 0030
C                                                                       RKF 0040
C     WRITTEN BY H.A.WATTS AND L.F.SHAMPINE                             RKF 0050
C                   SANDIA LABORATORIES                                 RKF 0060
C                  ALBUQUERQUE,NEW MEXICO                               RKF 0070
C                                                                       RKF 0080
C    RKF45 IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND MILDLY STIFF    RKF 0090
C    DIFFERENTIAL EQUATIONS WHEN DERIVATIVE EVALUATIONS ARE INEXPENSIVE.RKF 0100
C    RKF45 SHOULD GENERALLY NOT BE USED WHEN THE USER IS DEMANDING      RKF 0110
C    HIGH ACCURACY.                                                     RKF 0120
C                                                                       RKF 0130
C***********************************************************************RKF 0140
C ABSTRACT                                                              RKF 0150
C***********************************************************************RKF 0160
C                                                                       RKF 0170
C    SUBROUTINE  RKF45  INTEGRATES A SYSTEM OF NEQN FIRST ORDER         RKF 0180
C    ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM                        RKF 0190
C             DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN))                     RKF 0200
C              WHERE THE Y(I) ARE GIVEN AT T .                          RKF 0210
C    TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TOUT BUT ITRKF 0220
C    CAN BE USED AS A ONE-STEP INTEGRATOR TO ADVANCE THE SOLUTION A     RKF 0230
C    SINGLE STEP IN THE DIRECTION OF TOUT.  ON RETURN THE PARAMETERS IN RKF 0240
C    THE CALL LIST ARE SET FOR CONTINUING THE INTEGRATION. THE USER HAS RKF 0250
C    ONLY TO CALL RKF45 AGAIN (AND PERHAPS DEFINE A NEW VALUE FOR TOUT).RKF 0260
C    ACTUALLY, RKF45 IS AN INTERFACING ROUTINE WHICH CALLS SUBROUTINE   RKF 0270
C    RKFS FOR THE SOLUTION.  RKFS IN TURN CALLS SUBROUTINE  FEHL WHICH  RKF 0280
C    COMPUTES AN APPROXIMATE SOLUTION OVER ONE STEP.                    RKF 0290
C                                                                       RKF 0300
C    RKF45  USES THE RUNGE-KUTTA-FEHLBERG (4,5)  METHOD DESCRIBED       RKF 0310
C    IN THE REFERENCE                                                   RKF 0320
C    E.FEHLBERG , LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITH STEPSIZERKF 0330
C                 CONTROL , NASA TR R-315                               RKF 0340
C                                                                       RKF 0350
C    THE PERFORMANCE OF RKF45 IS ILLUSTRATED IN THE REFERENCE           RKF 0360
C    L.F.SHAMPINE,H.A.WATTS,S.DAVENPORT, SOLVING NON-STIFF ORDINARY     RKF 0370
C                 DIFFERENTIAL EQUATIONS-THE STATE OF THE ART ,         RKF 0380
C                 SIAM REVIEW 18(1976),PP.376-411                       RKF 0390
C                                                                       RKF 0400
C                                                                       RKF 0410
C    THE PARAMETERS REPRESENT-                                          RKF 0420
C      F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT RKF 0430
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     RKF 0440
C      Y(*) -- SOLUTION VECTOR AT T                                     RKF 0450
C      T -- INDEPENDENT VARIABLE                                        RKF 0460
C      TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED                RKF 0470
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCALRKF 0480
C            ERROR TEST. AT EACH STEP THE CODE REQUIRES THAT            RKF 0490
C                 ABS(LOCAL ERROR) .LE. RELERR*ABS(Y) + ABSERR          RKF 0500
C            FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS RKF 0510
C      IFLAG -- INDICATOR FOR STATUS OF INTEGRATION                     RKF 0520
C      WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO RKF45 WHICH IS  RKF 0530
C            NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED        RKF 0540
C            AT LEAST  3+6*NEQN                                         RKF 0550
C      IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL TO   RKF 0560
C            RKF45 WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE     RKF 0570
C            DIMENSIONED AT LEAST  5                                    RKF 0580
C                                                                       RKF 0590
C                                                                       RKF 0600
C***********************************************************************RKF 0610
C  FIRST CALL TO RKF45                                                  RKF 0620
C***********************************************************************RKF 0630
C                                                                       RKF 0640
C    THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYSRKF 0650
C    IN THE CALL LIST  -      Y(NEQN) , WORK(3+6*NEQN) , IWORK(5)  ,    RKF 0660
C    DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE F(T,Y,YP) ANDRKF 0670
C    INITIALIZE THE FOLLOWING PARAMETERS-                               RKF 0680
C                                                                       RKF 0690
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED.  (NEQN .GE. 1)     RKF 0700
C      Y(*) -- VECTOR OF INITIAL CONDITIONS                             RKF 0710
C      T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE          RKF 0720
C      TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED.               RKF 0730
C            T=TOUT IS ALLOWED ON THE FIRST CALL ONLY, IN WHICH CASE    RKF 0740
C            RKF45 RETURNS WITH IFLAG=2 IF CONTINUATION IS POSSIBLE.    RKF 0750
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES    RKF 0760
C            WHICH MUST BE NON-NEGATIVE. RELERR MUST BE A VARIABLE WHILERKF 0770
C            ABSERR MAY BE A CONSTANT. THE CODE SHOULD NORMALLY NOT BE  RKF 0780
C            USED WITH RELATIVE ERROR CONTROL SMALLER THAN ABOUT 1.E-8 .RKF 0790
C            TO AVOID LIMITING PRECISION DIFFICULTIES THE CODE REQUIRES RKF 0800
C            RELERR TO BE LARGER THAN AN INTERNALLY COMPUTED RELATIVE   RKF 0810
C            ERROR PARAMETER WHICH IS MACHINE DEPENDENT. IN PARTICULAR, RKF 0820
C            PURE ABSOLUTE ERROR IS NOT PERMITTED. IF A SMALLER THAN    RKF 0830
C            ALLOWABLE VALUE OF RELERR IS ATTEMPTED, RKF45 INCREASES    RKF 0840
C            RELERR APPROPRIATELY AND RETURNS CONTROL TO THE USER BEFORERKF 0850
C            CONTINUING THE INTEGRATION.                                RKF 0860
C      IFLAG -- +1,-1  INDICATOR TO INITIALIZE THE CODE FOR EACH NEW    RKF 0870
C            PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG=-1  RKF 0880
C            ONLY WHEN ONE-STEP INTEGRATOR CONTROL IS ESSENTIAL. IN THISRKF 0890
C            CASE, RKF45 ATTEMPTS TO ADVANCE THE SOLUTION A SINGLE STEP RKF 0900
C            IN THE DIRECTION OF TOUT EACH TIME IT IS CALLED. SINCE THISRKF 0910
C            MODE OF OPERATION RESULTS IN EXTRA COMPUTING OVERHEAD, IT  RKF 0920
C            SHOULD BE AVOIDED UNLESS NEEDED.                           RKF 0930
C                                                                       RKF 0940
C                                                                       RKF 0950
C***********************************************************************RKF 0960
C  OUTPUT FROM RKF45                                                    RKF 0970
C***********************************************************************RKF 0980
C                                                                       RKF 0990
C      Y(*) -- SOLUTION AT T                                            RKF 1000
C      T -- LAST POINT REACHED IN INTEGRATION.                          RKF 1010
C      IFLAG = 2 -- INTEGRATION REACHED TOUT.INDICATES SUCCESSFUL RETURNRKF 1020
C                   AND IS THE NORMAL MODE FOR CONTINUING INTEGRATION.  RKF 1030
C            =-2 -- A SINGLE SUCCESSFUL STEP IN THE DIRECTION OF TOUT   RKF 1040
C                   HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING          RKF 1050
C                   INTEGRATION ONE STEP AT A TIME.                     RKF 1060
C            = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE RELATIVE ERRORRKF 1070
C                   TOLERANCE WAS TOO SMALL. RELERR HAS BEEN INCREASED  RKF 1080
C                   APPROPRIATELY FOR CONTINUING.                       RKF 1090
C            = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN     RKF 1100
C                   3000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS       RKF 1110
C                   IS APPROXIMATELY 500 STEPS.                         RKF 1120
C            = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION      RKF 1130
C                   VANISHED MAKING A PURE RELATIVE ERROR TEST          RKF 1140
C                   IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE.   RKF 1150
C                   USING THE ONE-STEP INTEGRATION MODE FOR ONE STEP    RKF 1160
C                   IS A GOOD WAY TO PROCEED.                           RKF 1170
C            = 6 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED     RKF 1180
C                   ACCURACY COULD NOT BE ACHIEVED USING SMALLEST       RKF 1190
C                   ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR    RKF 1200
C                   TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE       RKF 1210
C                   ATTEMPTED.                                          RKF 1220
C            = 7 -- IT IS LIKELY THAT RKF45 IS INEFFICIENT FOR SOLVING  RKF 1230
C                   THIS PROBLEM. TOO MUCH OUTPUT IS RESTRICTING THE    RKF 1240
C                   NATURAL STEPSIZE CHOICE. USE THE ONE-STEP INTEGRATORRKF 1250
C                   MODE.                                               RKF 1260
C            = 8 -- INVALID INPUT PARAMETERS                            RKF 1270
C                   THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS    RKF 1280
C                   SATISFIED -   NEQN .LE. 0                           RKF 1290
C                                 T=TOUT  AND  IFLAG .NE. +1 OR -1      RKF 1300
C                                 RELERR OR ABSERR .LT. 0.              RKF 1310
C                                 IFLAG .EQ. 0  OR  .LT. -2  OR  .GT. 8 RKF 1320
C      WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO INTEREST  RKF 1330
C                   TO THE USER BUT NECESSARY FOR SUBSEQUENT CALLS.     RKF 1340
C                   WORK(1),...,WORK(NEQN) CONTAIN THE FIRST DERIVATIVESRKF 1350
C                   OF THE SOLUTION VECTOR Y AT T. WORK(NEQN+1) CONTAINSRKF 1360
C                   THE STEPSIZE H TO BE ATTEMPTED ON THE NEXT STEP.    RKF 1370
C                   IWORK(1) CONTAINS THE DERIVATIVE EVALUATION COUNTER.RKF 1380
C                                                                       RKF 1390
C                                                                       RKF 1400
C***********************************************************************RKF 1410
C  SUBSEQUENT CALLS TO RKF45                                            RKF 1420
C***********************************************************************RKF 1430
C                                                                       RKF 1440
C    SUBROUTINE RKF45 RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE   RKF 1450
C    THE INTEGRATION. IF THE INTEGRATION REACHED TOUT,THE USER NEED ONLYRKF 1460
C    DEFINE A NEW TOUT AND CALL RKF45 AGAIN. IN THE ONE-STEP INTEGRATOR RKF 1470
C    MODE (IFLAG=-2) THE USER MUST KEEP IN MIND THAT EACH STEP TAKEN IS RKF 1480
C    IN THE DIRECTION OF THE CURRENT TOUT. UPON REACHING TOUT (INDICATEDRKF 1490
C    BY CHANGING IFLAG TO 2),THE USER MUST THEN DEFINE A NEW TOUT AND   RKF 1500
C    RESET IFLAG TO -2 TO CONTINUE IN THE ONE-STEP INTEGRATOR MODE.     RKF 1510
C                                                                       RKF 1520
C    IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS TO   RKF 1530
C    CONTINUE (IFLAG=3,4 CASES), HE JUST CALLS RKF45 AGAIN. WITH IFLAG=3RKF 1540
C    THE RELERR PARAMETER HAS BEEN ADJUSTED APPROPRIATELY FOR CONTINUINGRKF 1550
C    THE INTEGRATION. IN THE CASE OF IFLAG=4 THE FUNCTION COUNTER WILL  RKF 1560
C    BE RESET TO 0 AND ANOTHER 3000 FUNCTION EVALUATIONS ARE ALLOWED.   RKF 1570
C                                                                       RKF 1580
C    HOWEVER,IN THE CASE IFLAG=5, THE USER MUST FIRST ALTER THE ERROR   RKF 1590
C    CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE INTEGRATION CAN RKF 1600
C    PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED.                   RKF 1610
C                                                                       RKF 1620
C    ALSO,IN THE CASE IFLAG=6, IT IS NECESSARY FOR THE USER TO RESET    RKF 1630
C    IFLAG TO 2 (OR -2 WHEN THE ONE-STEP INTEGRATION MODE IS BEING USED)RKF 1640
C    AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH BEFORE THE      RKF 1650
C    INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, EXECUTION WILL  RKF 1660
C    BE TERMINATED. THE OCCURRENCE OF IFLAG=6 INDICATES A TROUBLE SPOT  RKF 1670
C    %SOLUTION IS CHANGING RAPIDLY,SINGULARITY MAY BE PRESENT) AND IT   RKF 1680
C    OFTEN IS INADVISABLE TO CONTINUE.                                  RKF 1690
C                                                                       RKF 1700
C    IF IFLAG=7 IS ENCOUNTERED, THE USER SHOULD USE THE ONE-STEP        RKF 1710
C    INTEGRATION MODE WITH THE STEPSIZE DETERMINED BY THE CODE OR       RKF 1720
C    CONSIDER SWITCHING TO THE ADAMS CODES DE/STEP,INTRP. IF THE USER   RKF 1730
C    INSISTS UPON CONTINUING THE INTEGRATION WITH RKF45, HE MUST RESET  RKF 1740
C    IFLAG TO 2 BEFORE CALLING RKF45 AGAIN. OTHERWISE,EXECUTION WILL BE RKF 1750
C    TERMINATED.                                                        RKF 1760
C                                                                       RKF 1770
C    IF IFLAG=8 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS    RKF 1780
C    THE INVALID INPUT PARAMETERS ARE CORRECTED.                        RKF 1790
C                                                                       RKF 1800
C    IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN INFORMATION  RKF 1810
C    REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, WORK AND IWORK   RKF 1820
C    SHOULD NOT BE ALTERED.                                             RKF 1830
C                                                                       RKF 1840
C***********************************************************************RKF 1850
C                                                                       RKF 1860
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                RKF 1865
      DIMENSION Y(NEQN),WORK(1),IWORK(5)                                RKF 1870
C                                                                       RKF 1880
      EXTERNAL F                                                        RKF 1890
C                                                                       RKF 1900
C                                                                       RKF 1910
C     COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY               RKF 1920
C                                                                       RKF 1930
      K1M=NEQN+1                                                        RKF 1940
      K1=K1M+1                                                          RKF 1950
      K2=K1+NEQN                                                        RKF 1960
      K3=K2+NEQN                                                        RKF 1970
      K4=K3+NEQN                                                        RKF 1980
      K5=K4+NEQN                                                        RKF 1990
      K6=K5+NEQN                                                        RKF 2000
C                                                                       RKF 2010
C***********************************************************************RKF 2020
C     THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG       RKF 2030
C     CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE       RKF 2040
C     ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER,        RKF 2050
C     HE MUST USE RKFS DIRECTLY.                                        RKF 2060
C***********************************************************************RKF 2070
C                                                                       RKF 2080
      CALL RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK(1),WORK(K1M),  RKF 2090
     1          WORK(K1),WORK(K2),WORK(K3),WORK(K4),WORK(K5),WORK(K6),  RKF 2100
     2          WORK(K6+1),IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5),RKF 2110
     3          R,D)
C                                                                       RKF 2120
      RETURN                                                            RKF 2130
      END                                                               RKF 2140

      SUBROUTINE RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,YP,H,F1,F2,F3,RKFS0010
     1                F4,F5,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,KFLAG,R,D)   RKFS0020
C                                                                       RKFS0030
C     FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD                    RKFS0040
C                                                                       RKFS0050
C***********************************************************************RKFS0060
C                                                                       RKFS0070
C     RKFS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL     RKFS0080
C     EQUATIONS AS DESCRIBED IN THE COMMENTS FOR RKF45 .                RKFS0090
C     THE ARRAYS YP,F1,F2,F3,F4,AND F5 (OF DIMENSION AT LEAST NEQN) AND RKFS0100
C     THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE USED RKFS0110
C     INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO ELIMINATE   RKFS0120
C     LOCAL RETENTION OF VARIABLES BETWEEN CALLS. ACCORDINGLY, THEY     RKFS0130
C     SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE INTEREST ARE             RKFS0140
C         YP - DERIVATIVE OF SOLUTION VECTOR AT T                       RKFS0150
C         H  - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP     RKFS0160
C         NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION EVALUATIONS RKFS0170
C                                                                       RKFS0180
C***********************************************************************RKFS0190
C                                                                       RKFS0200
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                RKFS0205
      LOGICAL HFAILD,OUTPUT                                             RKFS0210
C                                                                       RKFS0220
      DIMENSION Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),   RKFS0230
     1          F5(NEQN)                                                RKFS0240
      PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,C26=26.0D0,CP2=0.2D0,    RKFS0241
     1   CP5=0.5D0,C2090=2090.D0,C21970=21970.D0,C15048=15048.D0,       RKFS0242
     2   C22528=22528.D0,C27360=27360.D0,C7524C=752400.D0,CP1=0.1D0,    RKFS0243
     3   CP9=0.9D0,C59049=59049.D0,FIVE=5.0D0,C1P889=1.889568D-4)       RKFS0244
C                                                                       RKFS0250
      EXTERNAL F                                                        RKFS0260
C                                                                       RKFS0270
C***********************************************************************RKFS0280
C                                                                       RKFS0290
C  THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE VALUE    RKFS0300
C  REPRESENTABLE IN THE MACHINE SUCH THAT  1.+ U .GT. 1.                RKFS0310
C                  VALUES TO BE USED ARE                                RKFS0320
C              U = 9.5E-7          FOR IBM 360/370                      RKFS0330
C              U = 1.5E-8          FOR UNIVAC 1108                      RKFS0340
C              U = 7.5E-9          FOR PDP-10                           RKFS0350
C              U = 7.1E-15         FOR CDC 6000  SERIES                 RKFS0360
C              U = 2.2E-16         FOR IBM 360/370  DOUBLE PRECISION    RKFS0370
C              U = 1.4E-14         FOR PRIME 650/750 DOUBLE PRECISION   RKFS0375
C                                                                       RKFS0380
      DATA U/1.4D-14/                                                   RKFS0390
C                                                                       RKFS0400
C***********************************************************************RKFS0410
C                                                                       RKFS0420
C  REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE       RKFS0430
C  INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL         RKFS0440
C  GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING      RKFS0450
C  PRECISION ON COMPUTERS WITH LONG WORDLENGTHS.                        RKFS0460
C                                                                       RKFS0470
      DATA REMIN/1.D-12/                                                RKFS0480
C                                                                       RKFS0490
C***********************************************************************RKFS0500
C                                                                       RKFS0510
C     THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER               RKFS0520
C     OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.               RKFS0530
C     AS SET,THIS CORRESPONDS TO ABOUT 500 STEPS.                       RKFS0540
C                                                                       RKFS0550
      DATA MAXNFE/3000/                                                 RKFS0560
C                                                                       RKFS0570
C***********************************************************************RKFS0580
C                                                                       RKFS0590
C                                                                       RKFS0600
C     CHECK INPUT PARAMETERS                                            RKFS0610
C                                                                       RKFS0620
C                                                                       RKFS0630
      IF (NEQN .LT. 1) GO TO 10                                         RKFS0640
      IF ((RELERR .LT. ZERO)  .OR.  (ABSERR .LT. ZERO)) GO TO 10        RKFS0650
      MFLAG=IABS(IFLAG)                                                 RKFS0660
      IF ((MFLAG .GE. 1)  .AND.  (MFLAG .LE. 8)) GO TO 20               RKFS0670
C                                                                       RKFS0680
C     INVALID INPUT                                                     RKFS0690
   10 IFLAG=8                                                           RKFS0700
      RETURN                                                            RKFS0710
C                                                                       RKFS0720
C     IS THIS THE FIRST CALL                                            RKFS0730
   20 IF (MFLAG .EQ. 1) GO TO 50                                        RKFS0740
C                                                                       RKFS0750
C     CHECK CONTINUATION POSSIBILITIES                                  RKFS0760
C                                                                       RKFS0770
      IF ((T .EQ. TOUT) .AND. (KFLAG .NE. 3)) GO TO 10                  RKFS0780
      IF (MFLAG .NE. 2) GO TO 25                                        RKFS0790
C                                                                       RKFS0800
C     IFLAG = +2 OR -2                                                  RKFS0810
      IF (KFLAG .EQ. 3) GO TO 45                                        RKFS0820
      IF (INIT .EQ. 0) GO TO 45                                         RKFS0830
      IF (KFLAG .EQ. 4) GO TO 40                                        RKFS0840
      IF ((KFLAG .EQ. 5)  .AND.  (ABSERR .EQ. ZERO)) GO TO 30           RKFS0850
      IF ((KFLAG .EQ. 6)  .AND.  (RELERR .LE. SAVRE)  .AND.             RKFS0860
     1    (ABSERR .LE. SAVAE)) GO TO 30                                 RKFS0870
      GO TO 50                                                          RKFS0880
C                                                                       RKFS0890
C     IFLAG = 3,4,5,6,7 OR 8                                            RKFS0900
   25 IF (IFLAG .EQ. 3) GO TO 45                                        RKFS0910
      IF (IFLAG .EQ. 4) GO TO 40                                        RKFS0920
      IF ((IFLAG .EQ. 5) .AND. (ABSERR .GT. ZERO)) GO TO 45             RKFS0930
C                                                                       RKFS0940
C     INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO     RKFS0950
C     THE INSTRUCTIONS PERTAINING TO IFLAG=5,6,7 OR 8                   RKFS0960
   30 STOP                                                              RKFS0970
C                                                                       RKFS0980
C***********************************************************************RKFS0990
C                                                                       RKFS1000
C     RESET FUNCTION EVALUATION COUNTER                                 RKFS1010
   40 NFE=0                                                             RKFS1020
      IF (MFLAG .EQ. 2) GO TO 50                                        RKFS1030
C                                                                       RKFS1040
C     RESET FLAG VALUE FROM PREVIOUS CALL                               RKFS1050
   45 IFLAG=JFLAG                                                       RKFS1060
      IF (KFLAG .EQ. 3) MFLAG=IABS(IFLAG)                               RKFS1070
C                                                                       RKFS1080
C     SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT   RKFS1090
C     INPUT CHECKING                                                    RKFS1100
   50 JFLAG=IFLAG                                                       RKFS1110
      KFLAG=0                                                           RKFS1120
C                                                                       RKFS1130
C     SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS     RKFS1140
      SAVRE=RELERR                                                      RKFS1150
      SAVAE=ABSERR                                                      RKFS1160
C                                                                       RKFS1170
C     RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS      RKFS1180
C     2U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING FROM    RKFS1190
C     IMPOSSIBLE ACCURACY REQUESTS                                      RKFS1200
C                                                                       RKFS1210
      RER=TWO*U+REMIN                                                   RKFS1220
      IF (RELERR .GE. RER) GO TO 55                                     RKFS1230
C                                                                       RKFS1240
C     RELATIVE ERROR TOLERANCE TOO SMALL                                RKFS1250
      RELERR=RER                                                        RKFS1260
      IFLAG=3                                                           RKFS1270
      KFLAG=3                                                           RKFS1280
      RETURN                                                            RKFS1290
C                                                                       RKFS1300
   55 U26=C26*U                                                         RKFS1310
C                                                                       RKFS1320
      DT=TOUT-T                                                         RKFS1330
C                                                                       RKFS1340
      IF (MFLAG .EQ. 1) GO TO 60                                        RKFS1350
      IF (INIT .EQ. 0) GO TO 65                                         RKFS1360
      GO TO 80                                                          RKFS1370
C                                                                       RKFS1380
C                                                                       RKFS1390
C***********************************************************************RKFS1400
C                                                                       RKFS1410
C     INITIALIZATION --                                                 RKFS1420
C                       SET INITIALIZATION COMPLETION INDICATOR,INIT    RKFS1430
C                       SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP    RKFS1440
C                       EVALUATE INITIAL DERIVATIVES                    RKFS1450
C                       SET COUNTER FOR FUNCTION EVALUATIONS,NFE        RKFS1460
C                       ESTIMATE STARTING STEPSIZE                      RKFS1470
C                                                                       RKFS1480
   60 INIT=0                                                            RKFS1490
      KOP=0                                                             RKFS1500
C                                                                       RKFS1510
      A=T                                                               RKFS1520
      CALL F(A,Y,YP,R,D)                                                RKFS1530
      NFE=1                                                             RKFS1540
      IF (T .NE. TOUT) GO TO 65                                         RKFS1550
      IFLAG=2                                                           RKFS1560
      RETURN                                                            RKFS1570
C                                                                       RKFS1580
C                                                                       RKFS1590
   65 INIT=1                                                            RKFS1600
      H=DABS(DT)                                                        RKFS1610
      TOLN=ZERO                                                         RKFS1620
      DO 70 K=1,NEQN                                                    RKFS1630
        TOL=RELERR*DABS(Y(K))+ABSERR                                    RKFS1640
        IF (TOL .LE. ZERO) GO TO 70                                     RKFS1650
        TOLN=TOL                                                        RKFS1660
        YPK=DABS(YP(K))                                                 RKFS1670
        IF (YPK*H**5 .GT. TOL) H=(TOL/YPK)**CP2                         RKFS1680
   70 CONTINUE                                                          RKFS1690
      IF (TOLN .LE. ZERO) H=ZERO                                        RKFS1700
      H=DMAX1(H,U26*DMAX1(DABS(T),DABS(DT)))                            RKFS1710
      JFLAG=ISIGN(2,IFLAG)                                              RKFS1720
C                                                                       RKFS1730
C                                                                       RKFS1740
C***********************************************************************RKFS1750
C                                                                       RKFS1760
C     SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT      RKFS1770
C                                                                       RKFS1780
   80 H=DSIGN(H,DT)                                                     RKFS1790
C                                                                       RKFS1800
C     TEST TO SEE IF RKF45 IS BEING SEVERELY IMPACTED BY TOO MANY       RKFS1810
C     OUTPUT POINTS                                                     RKFS1820
C                                                                       RKFS1830
      IF (DABS(H) .GE. TWO*DABS(DT)) KOP=KOP+1                          RKFS1840
      IF (KOP .NE. 100) GO TO 85                                        RKFS1850
C                                                                       RKFS1860
C     UNNECESSARY FREQUENCY OF OUTPUT                                   RKFS1870
      KOP=0                                                             RKFS1880
      IFLAG=7                                                           RKFS1890
      RETURN                                                            RKFS1900
C                                                                       RKFS1910
   85 IF (DABS(DT) .GT. U26*DABS(T)) GO TO 95                           RKFS1920
C                                                                       RKFS1930
C     IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN               RKFS1940
C                                                                       RKFS1950
      DO 90 K=1,NEQN                                                    RKFS1960
   90   Y(K)=Y(K)+DT*YP(K)                                              RKFS1970
      A=TOUT                                                            RKFS1980
      CALL F(A,Y,YP,R,D)                                                RKFS1990
      NFE=NFE+1                                                         RKFS2000
      GO TO 300                                                         RKFS2010
C                                                                       RKFS2020
C                                                                       RKFS2030
C     INITIALIZE OUTPUT POINT INDICATOR                                 RKFS2040
C                                                                       RKFS2050
   95 OUTPUT= .FALSE.                                                   RKFS2060
C                                                                       RKFS2070
C     TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION,     RKFS2080
C     SCALE THE ERROR TOLERANCES                                        RKFS2090
C                                                                       RKFS2100
      SCALE=TWO/RELERR                                                  RKFS2110
      AE=SCALE*ABSERR                                                   RKFS2120
C                                                                       RKFS2130
C                                                                       RKFS2140
C***********************************************************************RKFS2150
C***********************************************************************RKFS2160
C     STEP BY STEP INTEGRATION                                          RKFS2170
C                                                                       RKFS2180
  100 HFAILD= .FALSE.                                                   RKFS2190
      NS=NS+1                                                           RKFS2200
C                                                                       RKFS2210
C     SET SMALLEST ALLOWABLE STEPSIZE                                   RKFS2220
C                                                                       RKFS2230
      HMIN=U26*DABS(T)                                                  RKFS2240
C                                                                       RKFS2250
C     ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT.             RKFS2260
C     LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE AND RKFS2270
C     THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE.              RKFS2280
C                                                                       RKFS2290
      DT=TOUT-T                                                         RKFS2300
      IF (DABS(DT) .GE. TWO*DABS(H)) GO TO 200                          RKFS2310
      IF (DABS(DT) .GT. DABS(H)) GO TO 150                              RKFS2320
C                                                                       RKFS2330
C     THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE     RKFS2340
C     OUTPUT POINT                                                      RKFS2350
C                                                                       RKFS2360
      OUTPUT= .TRUE.                                                    RKFS2370
      H=DT                                                              RKFS2380
      GO TO 200                                                         RKFS2390
C                                                                       RKFS2400
  150 H=CP5*DT                                                          RKFS2410
C                                                                       RKFS2420
C                                                                       RKFS2430
C                                                                       RKFS2440
C***********************************************************************RKFS2450
C     CORE INTEGRATOR FOR TAKING A SINGLE STEP                          RKFS2460
C***********************************************************************RKFS2470
C     THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW IN   RKFS2480
C     COMPUTING THE ERROR TOLERANCE FUNCTION ET.                        RKFS2490
C     TO AVOID PROBLEMS WITH ZERO CROSSINGS,RELATIVE ERROR IS MEASURED  RKFS2500
C     USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION AT THE        RKFS2510
C     BEGINNING AND END OF A STEP.                                      RKFS2520
C     THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF    RKFS2530
C     SIGNIFICANCE.                                                     RKFS2540
C     TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED          RKFS2550
C     TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T.                 RKFS2560
C     PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO    RKFS2570
C     SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE      RKFS2580
C     CHATTERING ON PROBLEMS HAVING DISCONTINUITIES.                    RKFS2590
C     TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE STEPSIZE  RKFS2600
C     IT ESTIMATES WILL SUCCEED.                                        RKFS2610
C     AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE FOR RKFS2620
C     THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE EFFICIENT ON    RKFS2630
C     PROBLEMS HAVING DISCONTINUITIES AND MORE EFFECTIVE IN GENERAL     RKFS2640
C     SINCE LOCAL EXTRAPOLATION IS BEING USED AND EXTRA CAUTION SEEMS   RKFS2650
C     WARRANTED.                                                        RKFS2660
C***********************************************************************RKFS2670
C                                                                       RKFS2680
C                                                                       RKFS2690
C     TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS.                   RKFS2700
C     IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H              RKFS2710
C                                                                       RKFS2720
  200 IF (NFE .LE. MAXNFE) GO TO 220                                    RKFS2730
C                                                                       RKFS2740
C     TOO MUCH WORK                                                     RKFS2750
      IFLAG=4                                                           RKFS2760
      KFLAG=4                                                           RKFS2770
      RETURN                                                            RKFS2780
C                                                                       RKFS2790
C     ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H         RKFS2800
C                                                                       RKFS2810
  220 CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1,R,D)                  RKFS2820
      NFE=NFE+5                                                         RKFS2830
C                                                                       RKFS2840
C     COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR ESTIMATESRKFS2850
C     AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE ERROR IS     RKFS2860
C     MEASURED WITH RESPECT TO THE AVERAGE OF THE MAGNITUDES OF THE     RKFS2870
C     SOLUTION AT THE BEGINNING AND END OF THE STEP.                    RKFS2880
C                                                                       RKFS2890
      EEOET=ZERO                                                        RKFS2900
      DO 250 K=1,NEQN                                                   RKFS2910
        ET=DABS(Y(K))+DABS(F1(K))+AE                                    RKFS2920
        IF (ET .GT. ZERO) GO TO 240                                     RKFS2930
C                                                                       RKFS2940
C       INAPPROPRIATE ERROR TOLERANCE                                   RKFS2950
        IFLAG=5                                                         RKFS2960
        KFLAG=5                                                         RKFS2970
        RETURN                                                          RKFS2980
C                                                                       RKFS2990
  240   EE=DABS((-C2090*YP(K)+(C21970*F3(K)-C15048*F4(K)))+             RKFS3000
     1                        (C22528*F2(K)-C27360*F5(K)))              RKFS3010
  250   EEOET=DMAX1(EEOET,EE/ET)                                        RKFS3020
C                                                                       RKFS3030
      ESTTOL=DABS(H)*EEOET*SCALE/C7524C                                 RKFS3040
C                                                                       RKFS3050
      IF (ESTTOL .LE. ONE) GO TO 260                                    RKFS3060
C                                                                       RKFS3070
C                                                                       RKFS3080
C     UNSUCCESSFUL STEP                                                 RKFS3090
C                       REDUCE THE STEPSIZE , TRY AGAIN                 RKFS3100
C                       THE DECREASE IS LIMITED TO A FACTOR OF 1/10     RKFS3110
C                                                                       RKFS3120
      HFAILD= .TRUE.                                                    RKFS3130
      NFS=NFS+1                                                         RKFS3140
      OUTPUT= .FALSE.                                                   RKFS3150
      S=CP1                                                             RKFS3160
      IF (ESTTOL .LT. C59049) S=CP9/ESTTOL**CP2                         RKFS3170
      H=S*H                                                             RKFS3180
      IF (DABS(H) .GT. HMIN) GO TO 200                                  RKFS3190
C                                                                       RKFS3200
C     REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE       RKFS3210
      IFLAG=6                                                           RKFS3220
      KFLAG=6                                                           RKFS3230
      RETURN                                                            RKFS3240
C                                                                       RKFS3250
C                                                                       RKFS3260
C     SUCCESSFUL STEP                                                   RKFS3270
C                        STORE SOLUTION AT T+H                          RKFS3280
C                        AND EVALUATE DERIVATIVES THERE                 RKFS3290
C                                                                       RKFS3300
  260 T=T+H                                                             RKFS3310
      DO 270 K=1,NEQN                                                   RKFS3320
  270   Y(K)=F1(K)                                                      RKFS3330
      A=T                                                               RKFS3340
      CALL F(A,Y,YP,R,D)                                                RKFS3350
      NFE=NFE+1                                                         RKFS3360
C                                                                       RKFS3370
C                                                                       RKFS3380
C                       CHOOSE NEXT STEPSIZE                            RKFS3390
C                       THE INCREASE IS LIMITED TO A FACTOR OF 5        RKFS3400
C                       IF STEP FAILURE HAS JUST OCCURRED, NEXT         RKFS3410
C                          STEPSIZE IS NOT ALLOWED TO INCREASE          RKFS3420
C                                                                       RKFS3430
      S=FIVE                                                            RKFS3440
      IF (ESTTOL .GT. C1P889) S=CP9/ESTTOL**CP2                         RKFS3450
      IF (HFAILD) S=DMIN1(S,ONE)                                        RKFS3460
      H=DSIGN(DMAX1(S*DABS(H),HMIN),H)                                  RKFS3470
C                                                                       RKFS3480
C***********************************************************************RKFS3490
C     END OF CORE INTEGRATOR                                            RKFS3500
C***********************************************************************RKFS3510
C                                                                       RKFS3520
C                                                                       RKFS3530
C                                                                       RKFS3540
C     SHOULD WE TAKE ANOTHER STEP                                       RKFS3550
C                                                                       RKFS3560
      IF (OUTPUT) GO TO 300                                             RKFS3570
      IF (IFLAG .GT. 0) GO TO 100                                       RKFS3580
C                                                                       RKFS3590
C***********************************************************************RKFS3600
C***********************************************************************RKFS3610
C                                                                       RKFS3620
C                                                                       RKFS3630
C     INTEGRATION SUCCESSFULLY COMPLETED                                RKFS3640
C                                                                       RKFS3650
C     ONE-STEP MODE                                                     RKFS3660
      IFLAG=-2                                                          RKFS3670
      RETURN                                                            RKFS3680
C                                                                       RKFS3690
C     INTERVAL MODE                                                     RKFS3700
  300 T=TOUT                                                            RKFS3710
      IFLAG=2                                                           RKFS3720
      RETURN                                                            RKFS3730
C                                                                       RKFS3740
      END                                                               RKFS3750

      SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S,R,D)             FEHL0010
C                                                                       FEHL0020
C     FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD                    FEHL0030
C                                                                       FEHL0040
C***********************************************************************FEHL0050
C    FEHL INTEGRATES A SYSTEM OF NEQN FIRST ORDER                       FEHL0060
C    ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM                        FEHL0070
C             DY(I)/DT=F(T,Y(1),---,Y(NEQN))                            FEHL0080
C    WHERE THE INITIAL VALUES Y(I) AND THE INITIAL DERIVATIVES          FEHL0090
C    YP(I) ARE SPECIFIED AT THE STARTING POINT T. FEHL ADVANCES         FEHL0100
C    THE SOLUTION OVER THE FIXED STEP H AND RETURNS                     FEHL0110
C    THE FIFTH ORDER (SIXTH ORDER ACCURATE LOCALLY) SOLUTION            FEHL0120
C    APPROXIMATION AT T+H IN ARRAY S(I).                                FEHL0130
C    F1,---,F5 ARE ARRAYS OF DIMENSION NEQN WHICH ARE NEEDED            FEHL0140
C    FOR INTERNAL STORAGE.                                              FEHL0150
C    THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF SIGNIFICANCE.    FEHL0160
C    FEHL SHOULD BE CALLED WITH AN H NOT SMALLER THAN 13 UNITS OF       FEHL0170
C    ROUNDOFF IN T SO THAT THE VARIOUS INDEPENDENT ARGUMENTS CAN BE     FEHL0180
C    DISTINGUISHED.                                                     FEHL0190
C***********************************************************************FEHL0200
C                                                                       FEHL0210
C                                                                       FEHL0220
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                FEHL0225
      DIMENSION Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),   FEHL0230
     1          F5(NEQN),S(NEQN)                                        FEHL0240
      PARAMETER(TWO=2.D0,THREE=3.D0,FOUR=4.D0,EIGHT=8.D0,C32=32.D0)     FEHL0245     
      EXTERNAL F                                                         FEHL0250
      CH=H/FOUR                                                         FEHL0260
      DO 221 K=1,NEQN                                                   FEHL0270
  221   F5(K)=Y(K)+CH*YP(K)                                             FEHL0280
      CALL F(T+CH,F5,F1,R,D)                                            FEHL0290
C                                                                       FEHL0300
      CH=THREE*H/C32                                                    FEHL0310
      DO 222 K=1,NEQN                                                   FEHL0320
  222   F5(K)=Y(K)+CH*(YP(K)+THREE*F1(K))                               FEHL0330
      CALL F(T+THREE*H/EIGHT,F5,F2,R,D)                                 FEHL0340
C                                                                       FEHL0350
      CH=H/2197.D0                                                      FEHL0360
      DO 223 K=1,NEQN                                                   FEHL0370
  223   F5(K)=Y(K)+CH*(1932.D0*YP(K)+(7296.D0*F2(K)-7200.D0*F1(K)))     FEHL0380
      CALL F(T+12.D0*H/13.D0,F5,F3,R,D)                                 FEHL0390
C                                                                       FEHL0400
      CH=H/4104.D0                                                      FEHL0410
      DO 224 K=1,NEQN                                                   FEHL0420
  224   F5(K)=Y(K)+CH*((8341.D0*YP(K)-845.D0*F3(K))+                    FEHL0430
     1                            (29440.D0*F2(K)-32832.D0*F1(K)))      FEHL0440
      CALL F(T+H,F5,F4,R,D)                                             FEHL0450
C                                                                       FEHL0460
      CH=H/20520.D0                                                     FEHL0470
      DO 225 K=1,NEQN                                                   FEHL0480
  225   F1(K)=Y(K)+CH*((-6080.D0*YP(K)+(9295.D0*F3(K)-5643.D0*F4(K)))+  FEHL0490
     1                             (41040.D0*F1(K)-28352.D0*F2(K)))     FEHL0500
      CALL F(T+H/TWO,F1,F5,R,D)                                         FEHL0510
C                                                                       FEHL0520
C     COMPUTE APPROXIMATE SOLUTION AT T+H                               FEHL0530
C                                                                       FEHL0540
      CH=H/7618050.D0                                                   FEHL0550
      DO 230 K=1,NEQN                                                   FEHL0560
  230   S(K)=Y(K)+CH*((902880.D0*YP(K)+(3855735.D0*F3(K)                FEHL0570
     1        -1371249.D0*F4(K)))+(3953664.D0*F2(K)+277020.D0*F5(K)))   FEHL0580
C                                                                       FEHL0590
      RETURN                                                            FEHL0600
      END                                                               FEHL0610






