      subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs,nd)
!
!     fifth-order RK
!
      implicit real *8 (a-h,o-z)
      parameter (nmax=100, pgrow=-.20,pshrnk=-.25,fcor=1.d0/15.,one=1., safety=.9, errcon=6.e-4,nd2=100)
      external derivs
      dimension y(nd),dydx(nd),yscal(nd),ytemp(nd2),ysav(nd2),dysav(nd2)
      xsav=x
      do 11 i=1,n
         ysav(i)=y(i)
         dysav(i)=dydx(i)
 11   continue
      h=htry
 1    hh=0.5*h
      call rk4(ysav,dysav,n,xsav,hh,ytemp,derivs,nd)
      x=xsav+hh
      call derivs(n,x,ytemp,dydx,nd)
      call rk4(ytemp,dydx,n,x,hh,y,derivs,nd)
      x=xsav+h
      if (x.eq.xsav) then
         write(*,*) 'stepsize not significant in rkqc',x
         stop
      endif
      call rk4(ysav,dysav,n,xsav,h,ytemp,derivs,nd)
      errmax=0.
      do 12 i=1,n
         ytemp(i)=y(i)-ytemp(i)
         dummy=abs(ytemp(i)/yscal(i))
         errmax=max(errmax,dummy)
 12   continue
      errmax=errmax/eps
      if(errmax.gt.one) then
         h=safety*h*(errmax**pshrnk)
         goto 1
      else
         hdid=h
         if (errmax.gt.errcon) then
            hnext=safety*h*(errmax**pgrow)
         else
            hnext=4.d0*h
         endif
      endif
      do 13 i=1,n
         y(i)=y(i)+ytemp(i)*fcor
 13   continue
      return
      end
!---------------------------------------------------------------------
      subroutine rk4(y,dydx,n,x,h,yout,derivs,nd)
!
!     rk4
!
      implicit real *8 (a-h,o-z)
      parameter (nmax=100,nd3=100)
      dimension y(nd),dydx(nd),yout(nd),yt(nd3),dyt(nd3),dym(nd3)
      external derivs
      hh=h*.5
      h6=h/6
      xh=x+hh
      do 11 i=1,n
         yt(i)=y(i)+hh*dydx(i)
 11   continue
      call derivs(n,xh,yt,dyt,nd)
      do 12 i=1,n
         yt(i)=y(i)+hh*dyt(i)
 12   continue
      call derivs(n,xh,yt,dym,nd)
      do 13 i=1,n
         yt(i)=y(i)+h*dym(i)
         dym(i)=dyt(i)+dym(i)
 13   continue
      call derivs(n,x+h,yt,dyt,nd)
      do 14 i=1,n
         yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2*dym(i))
 14   continue
      return
      end
