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
      IMPLICIT DOUBLE PRECISION(A-H,O-Z) 
      PARAMETER (ZERO=0.0D0)
      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='fehl45.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
      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
      WRITE(*,1000) NUMDE, XZERO, XEND, YZERO, H, IPRINT
1000  FORMAT(//,' EQUATION',I2,5X,'XZERO=',1PE9.2,5X,'B=',E9.2,
     *       5X,'YZERO=',E12.5,/,' STEPSIZE=',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
      WRITE(*,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
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BLOCKF/NUMDE
C
      GO TO(10,20,30), NUMDE
10    F = -Z
      RETURN
20    F = (Z + X*X - 2.0)/(X + 1.0)
      RETURN
30    F = -Z + 2.0D0*COS(X)
      RETURN
      END

      FUNCTION Y(X)
C     ----------
C
C     THIS GIVES THE TRUE SOLUTION OF
C     THE INITIAL VALUE PROBLEM.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BLOCKF/NUMDE
C
      GO TO(10,20,30), NUMDE
10    Y = EXP(-X)
      RETURN
20    X1 = X + 1.0
      Y = X*X - 2.0*(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 GIVEN (APPROXIMATE) VALUE OF THE SOLUTION
C         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
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D0)
      DIMENSION A(5),B(5,0:4),C(0:4),D(0:5),V(0:5)
C     THE FOLLOWING CONSTANTS ARE NEEDED IN EVALUATING
C     THE FEHLBERG FORMULA.
      PARAMETER(B30=1932.D0/2197.D0,B31=-7200.D0/2197.D0,
     *   B32=7296.D0/2197.D0,B40=439.D0/216.D0,B42=3680.D0/513.D0,
     *   B43=-845.D0/4104.D0,B50=-8.0D0/27.0D0,B52=-3544.D0/2565.D0,
     *   B53=1859.D0/4104.D0,B54=-11.D0/40.D0,A3=12.D0/13.D0,
     *   C0=25.D0/216.D0,C2=1408.D0/2565.D0,C3=2197.D0/4104.D0,
     *   D0=16.D0/135.D0,D2=6656.D0/12825.D0,D3=28561.D0/56430.D0,
     *   D5=2.D0/55.D0)
C     THE FOLLOWING ARRAYS ARE THE CONSTANTS IN THE FEHLBERG FORMULAS.
      DATA A/.25D0,.375D0,A3,1.0D0,0.5D0/,
     *     B/.25D0,.09375D0,B30,B40,B50,0.0D0,.28125D0,B31,-8.0D0,
     *       2.0D0,0.0D0,0.0D0,B32,B42,B52,0.0D0,0.0D0,0.0D0,B43,B53,
     *       0.0D0,0.0D0,0.0D0,0.0D0,B54/,
     *     C/C0,0.0D0,C2,C3,-.2D0/,
     *     D/D0,0.0D0,D2,D3,-.18D0,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