subroutine fehl(f,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s)
c
c     fehlberg fourth-fifth order runge-kutta method
c
c    fehl integrates a system of neqn first order
c    ordinary differential equations of the form
c             dy(i)/dt=f(t,y(1),---,y(neqn))
c    where the initial values y(i) and the initial derivatives
c    yp(i) are specified at the starting point t. fehl advances
c    the solution over the fixed step h and returns
c    the fifth order (sixth order accurate locally) solution
c    approximation at t+h in array s(i).
c    f1,---,f5 are arrays of dimension neqn which are needed
c    for internal storage.
c    the formulas have been grouped to control loss of significance.
c    fehl should be called with an h not smaller than 13 units of
c    roundoff in t so that the various independent arguments can be
c    distinguished.
c
c
      integer  neqn
      real  y(neqn),t,h,yp(neqn),f1(neqn),f2(neqn),
     1  f3(neqn),f4(neqn),f5(neqn),s(neqn)
c
      real  ch
      integer  k
c
      ch=h/4.0
      do 221 k=1,neqn
  221   f5(k)=y(k)+ch*yp(k)
      call f(t+ch,f5,f1)
c
      ch=3.0*h/32.0
      do 222 k=1,neqn
  222   f5(k)=y(k)+ch*(yp(k)+3.0*f1(k))
      call f(t+3.0*h/8.0,f5,f2)
c
      ch=h/2197.0
      do 223 k=1,neqn
  223   f5(k)=y(k)+ch*(1932.0*yp(k)+(7296.0*f2(k)-7200.0*f1(k)))
      call f(t+12.0*h/13.0,f5,f3)
c
      ch=h/4104.0
      do 224 k=1,neqn
  224   f5(k)=y(k)+ch*((8341.0*yp(k)-845.0*f3(k))+
     1                            (29440.0*f2(k)-32832.0*f1(k)))
      call f(t+h,f5,f4)
c
      ch=h/20520.0
      do 225 k=1,neqn
  225   f1(k)=y(k)+ch*((-6080.0*yp(k)+(9295.0*f3(k)-
     1         5643.0*f4(k)))+(41040.0*f1(k)-28352.0*f2(k)))
      call f(t+h/2.0,f1,f5)
c
c     compute approximate solution at t+h
c
      ch=h/7618050.0
      do 230 k=1,neqn
  230   s(k)=y(k)+ch*((902880.0*yp(k)+(3855735.0*f3(k)-
     1        1371249.0*f4(k)))+(3953664.0*f2(k)+
     2        277020.0*f5(k)))
c
      return
      end

.