subroutine fehl(f,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s)

      implicit double precision(a-h,o-z)
      real*8 y(8),yp(8),f1(8),f2(8),f3(8),f4(8),f5(8),s(8)

      ch=0.25*h
      do 221 k=1,neqn
        f5(k)=y(k)+ch*yp(k)
221   continue
      call f(t+0.25*h,f5,f1)
      ch=0.09375*h
      do 222 k=1,neqn
        f5(k)=y(k)+ch*(yp(k)+3.*f1(k))
222   continue
      call f(t+0.375*h,f5,f2) 
      ch=h/2197.
      do 223 k=1,neqn
        f5(k)=y(k)+ch*(1932.*yp(k)+(7296.*f2(k)-7200.*f1(k)))
223   continue
      call f(t+12./13.*h,f5,f3)
      ch=h/4104.
      do 224 k=1,neqn
        f5(k)=y(k)+ch*((8341.*yp(k)-845.*f3(k))+
     1                   (29440.*f2(k)-32832.*f1(k)))
224   continue
      call f(t+h,f5,f4)
      ch=h/20520.
      do 225 k=1,neqn
        f1(k)=y(k)+ch*((-6080.*yp(k)+(9295.*f3(k)-5643.*f4(k)))+
     1                             (41040.*f1(k)-28352.*f2(k)))
225   continue
      call f(t+0.5*h,f1,f5)
      ch=h/7618050. 
      do 230 k=1,neqn
        s(k)=y(k)+ch*((902880.*yp(k)+(3855735.*f3(k)-1371249.*f4(k)))+
     1                (3953664.*f2(k)+277020.*f5(k)))
230   continue

      return
      end 

************************************************************************