      program fehl45_test

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


