subroutine rkf(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)

      implicit double precision(a-h,o-z)

      real*8 y(8),work(51)
      real*8 yp(8),f1(8),f2(8),f3(8),f4(8),f5(8)
      integer iwork(5)
      external f

      k1m=neqn+1
      k1=k1m+1
      k2=k1+neqn
      k3=k2+neqn
      k4=k3+neqn
      k5=k4+neqn
      k6=k5+neqn
      do 10 i=1,neqn
        yp(i)=work(i)
10    continue
        h=work(k1m)
      do 20 j=1,neqn
        f1(j)=work(k1m+j)
20    continue
      do 30 k=1,neqn
        f2(k)=work(k2-1+k)
30    continue
      do 40 j=1,neqn
        f3(j)=work(k3-1+j)
40    continue
      do 50 k=1,neqn
        f4(k)=work(k4-1+k)
50    continue
      do 60 k=1,neqn
        f5(k)=work(k5-1+k)
60    continue
      savre=work(k6)
      savae=work(k6+1)
      call rkfs(f,neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1,f2,
     1 f3,f4,f5,savre,savae,iwork(1),iwork(2),iwork(3),iwork(4),
     2 iwork(5))
      work(k1m)=h
      work(k6)=savre
      work(k6+1)=savae
      do 70 k=1,neqn
        work(k)=yp(k)
70    continue
      do 80 j=1,neqn
        work(k1m+j)=f1(j)
80    continue
      do 90 k=1,neqn
        work(k2-1+k)=f2(k)
90    continue
      do 100 k=1,neqn
        work(k3-1+k)=f3(k)
100   continue
      do 110 k=1,neqn
        work(k4-1+k)=f4(k)
110   continue
      do 120 k=1,neqn
        work(k5-1+k)=f5(k)
120   continue
      return
      end

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