      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
c  Reference:
c
c    Kendall Atkinson, Weimin Han,
c    Elementary Numerical Analysis,
c    Wiley, 2004,
c    ISBN: 0471433373,
c    LC: QA297.A83.2004.
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. )
      parameter ( B31=-7200./2197. )
      parameter ( B32=7296./2197.,B40=439./216.,B42=3680./513.,
     *   B43=-845./4104.,B50=-8.0/27.0,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
