rkqs Subroutine

      SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)

      INTEGER n,NMAX
      double precision eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
      EXTERNAL derivs
      PARAMETER (NMAX=50)
!     CU USES derivs,rkck
      INTEGER i
      double precision errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,&
         PSHRNK,ERRCON
      PARAMETER (SAFETY=0.9d0,PGROW=-.2d0,PSHRNK=-.25d0,ERRCON=1.89d-4)

      h=htry

       if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x &
          .ne. x) then
         write(*,*) "RKQS BEFORE CALL RKCK: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         return
       end if

 1    call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
       if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x &
          .ne. x) then
         write(*,*) "RKQS AFTER CALL RKCK: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         return
       end if

      errmax=0.d0
      do 11 i=1,n
         errmax=max(errmax,abs(yerr(i)/yscal(i)))
 11   continue
      errmax=errmax/eps
      if(errmax.gt.1.d0)then
         h=SAFETY*h*(errmax**PSHRNK)
         if(h.lt.0.1d0*h)then
            h=.1d0*h
         end if
         xnew=x+h
         if(xnew.eq.x) call stop('stepsize underflow in rkqs')
         goto 1
      else
         if(errmax.gt.ERRCON)then
            hnext=SAFETY*h*(errmax**PGROW)
         else
            hnext=5.d0*h
         end if
         hdid=h
         x=x+h
         do 12 i=1,n
            y(i)=ytemp(i)
 12      continue
         return
      end if
      END subroutine rkqs