C     TITLE: DEMONSTRATION OF FEHL45.
C
C     THIS SOLVES THE INITIAL VALUE PROBLEM
C        Y'(X) = F(X,Y(X)) ,  X0 .LE. X .LE. B ,  Y(X0)=Y0.
C     THE FUNCTION F(X,Z) IS DEFINED BELOW, ALONG WITH THE TRUE
C     SOLUTION Y(X).  THE NUMBER OF THE PROBLEM TO BE SOLVED
C     IS SPECIFIED BY THE INPUT VARIABLE 'NUMDE', WHICH IS USED
C     IN THE FUNCTIONS 'F' AND 'Y'.  THE PROGRAM WILL REQUEST
C     THE PROBLEM PARAMETERS, ALONG WITH THE VALUES OF 'H' AND
C     'IPRINT'.  'H' IS THE STEPSIZE, AND 'IPRINT' IS THE NUMBER
C     OF STEPS BETWEEN EACH PRINTING OF THE SOLUTION.
C     USE H=0 AND NUMDE=0 TO STOP THE PROGRAM.
C
      PARAMETER(ZERO=0.0)
      COMMON/BLOCKF/NUMDE
C     *******************************************************
C     GIVE AN OUTPUT UNIT NUMBER, FOR A STORED OUTPUT FILE.
      DATA IOUT/8/
C     *******************************************************
      EXTERNAL F
C
      OPEN(IOUT, FILE='rkf45.output')
C     INPUT PROBLEM PARAMETERS.
10    PRINT *, ' WHICH DIFFERENTIAL EQUATION?'
      PRINT *, ' GIVE ZERO TO STOP.'
      READ *, NUMDE
      IF(NUMDE .EQ. 0) THEN
          CLOSE(IOUT)
          STOP
      END IF 
C
      PRINT *, ' GIVE DOMAIN [X0,B] OF SOLUTION.'
      READ *, XZERO, XEND
      PRINT *, ' WHAT IS Y0=Y(X0)?'
      READ *, YZERO
C
20    PRINT *, ' GIVE STEPSIZE H AND PRINT PARAMETER IPRINT.'
      PRINT *, ' LET H=0 TO TRY ANOTHER DIFFERENTIAL EQUATION.'
      READ *, H, IPRINT
      IF(H .EQ. ZERO) GO TO 10
C
C     INITIALIZE.
      X0 = XZERO
      Y0 = YZERO
      WRITE(IOUT,1000) NUMDE, XZERO, XEND, YZERO, H, IPRINT
1000  FORMAT(//,' EQUATION',I2,5X,'XZERO =',1PE9.2,5X,'B =',E9.2,
     *       5X,'YZERO =',E12.5,/,' STEPSIZ E=',E10.3,5X,
     *       'PRINT PARAMETER =',I3,/)
C
C     BEGIN MAIN LOOP FOR COMPUTING SOLUTION OF DIFFERENTIAL EQUATION.
30    DO 40 K=1,IPRINT
        X1 = X0  +H
        IF(X1 .GT. XEND) GO TO 20
        CALL FEHL45(F,Y0,X0,H,TRUN)
40      CONTINUE
C
C     CALCULATE ERROR AND PRINT RESULTS.
      TRUE = Y(X0)
      ERROR = TRUE - Y0
      WRITE(IOUT,1001) X0, Y0, ERROR, TRUN
1001  FORMAT(' X =',1PE10.3,5X,'Y(X) =',E17.10,5X,'ERROR =',E9.2,5X,
     *       'TRUN =',E9.2)
      GO TO 30
      END

      FUNCTION F(X,Z) 
C     ----------
C
C     THIS DEFINES THE RIGHT SIDE OF
C     THE DIFFERENTIAL EQUATION.
C
      PARAMETER (ONE=1.0, TWO=2.0)
      COMMON/BLOCKF/NUMDE
C
      GO TO(10,20,30), NUMDE
10    F = -Z
      RETURN
20    F = (Z + X*X - TWO)/(X + ONE)
      RETURN
30    F = -Z + TWO*COS(X)
      RETURN
      END

      FUNCTION Y(X)   
C     ----------
C
C     THIS GIVES THE TRUE SOLUTION OF
C     THE INITIAL VALUE PROBLEM.
C
      PARAMETER (ONE=1.0, TWO=2.0)
      COMMON/BLOCKF/NUMDE
C
      GO TO(10,20,30), NUMDE
10    Y = EXP(-X)
      RETURN
20    X1 = X + ONE
      Y = X*X - TWO*(X1*LOG(X1) - X1)
      RETURN
30    Y = SIN(X) + COS(X)
      RETURN
      END

      SUBROUTINE FEHL45(F,Y,X,H,TRUN)
C     -----------------
C
C     THIS COMPUTES A SINGLE STEP OF THE RUNGE-KUTTA-FEHLBERG METHOD
C     OF ORDER (4,5) FOR SOLVING THE DIFFERENTIAL EQUATION
C             Y'(X)=F(X,Y(X)).
C
C     INPUT:
C     (1) 'F' IS THE NAME OF A FUNCTION SUBPROGRAM FOR EVALUATING THE
C         DERIVATIVE F(X,Y) THAT DEFINES THE DIFFERENTIAL EQUATION.
C     (2) Y IS THE KNOWN VALUE OF THE SOLUTION AT THE POINT X.
C     (3) H IS THE STEPSIZE IN SOLVING THE DIFFERENTIAL EQUATION.
C         THE SOLUTION IS TO BE FOUND AT X+H.
C
C     OUTPUT:
C     (1) X WILL BE REPLACED BY X+H.
C     (2) Y WILL BE  REPLACED BY THE NUMERICAL SOLUTION BASED ON
C         USING THE FOURTH ORDER MEMBER OF THE FEHLBERG FORMULA
C         AT X+H.
C     (3) TRUN WILL BE SET TO THE ESTIMATED TRUNCATION ERROR IN THE
C         FOURTH ORDER SOLUTION Y; IT IS OBTAINED USING THE FIFTH
C         ORDER FEHLBERG SOLUTION.
C
C     TO SOLVE THE DIFFERENTIAL EQUATION ON A LONGER INTERVAL,
C     CALL THIS ROUTINE REPEATEDLY.
C

      DIMENSION A(5),B(5,0:4),C(0:4),D(0:5),V(0:5) 
      PARAMETER (ZERO=0.0)
C     THE FOLLOWING CONSTANTS ARE NEEDED IN EVALUATING
C     THE FEHLBERG FORMULA.
      PARAMETER(B30=1932./2197.,B31=-7200./2197.,
     *   B32=7296./2197.,B40=439./216.,B42=3680./513.,
     *   B43=-845./4104.,B50=-8.0D0/27.0D0,B52=-3544./2565.,
     *   B53=1859./4104.,B54=-11./40.,A3=12./13.,
     *   C0=25./216.,C2=1408./2565.,C3=2197./4104.,
     *   D0=16./135.,D2=6656./12825.,D3=28561./56430.,
     *   D5=2./55.)
C     THE FOLLOWING ARRAYS ARE THE CONSTANTS IN THE FEHLBERG FORMULAS.
      DATA A/.25,.375,A3,1.0,0.5/,
     *     B/.25,.09375,B30,B40,B50,0.0,.28125,B31,-8.0,
     *       2.0,0.0,0.0,B32,B42,B52,0.0,0.0,0.0,B43,B53,
     *       0.0,0.0,0.0,0.0,B54/,
     *     C/C0,0.0,C2,C3,-.2/,
     *     D/D0,0.0,D2,D3,-.18,D5/
C
C     EVALUATE THE DERIVATIVES.
      V(0) = F(X,Y)
      DO 20 I=1,5
        SUM  = ZERO
        DO 10 J=0,I-1
10        SUM = SUM + B(I,J)*V(J)
20      V(I) = F(X + A(I)*H, Y + H*SUM)
C
C     EVALUATE THE FOURTH AND FIFTH ORDER FORMULAS.
      FSUM = ZERO
      DO 30 I=0,4
30      FSUM = FSUM + C(I)*V(I)
      Y4 = Y + H*FSUM
      FSUM = ZERO
      DO 40 I=0,5
40      FSUM = FSUM + D(I)*V(I)
      Y5 = Y + H*FSUM
C
C     SET THE OUTPUT VARIABLES.
      X = X + H
      Y = Y4
      TRUN = Y5 - Y4
      RETURN
      END