      subroutine rkc4a(n,t0,tf,u0,nsteps,u)
C
C Subroutine rkc4a computes an ODE solution by a fixed step
C rkc4 method for a series of points along the solution by 
C repeatedly calling subroutine ssrkc4 for a single rkc4 step.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   tf       final value of independent variable
C
C   u0       initial condition vector of length n
C
C   nsteps   number of rkc4 steps
C
C   u        ODE solution vector of length n after
C            nsteps steps
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size arrays
      parameter(neq=250)
      dimension   e(neq)
      dimension u0(n), u(n)
C
C Integration step
      h=(tf-t0)/dfloat(nsteps)
C
C nsteps rkc4 steps
      do j=1,nsteps
C
C       Single rkc4 step
        call ssrkc4(n,t0,u0,h,t,u,e)
C
C       Reset base point values for the next rkc4 step
        do i=1,n
          u0(i)=u(i)
        end do
        t0=t
C
C     Next rkc4 step
      end do
C
C nsteps rkc4 steps completed
      return
C
C End of rkc4a
      end
      subroutine rkc4b(n,t0,tf,u0,nsteps,abserr,relerr,u)
C
C Subroutine rkc4b computes an ODE solution by a variable step
C classical fourth order RK method for a series of points along 
C the solution by repeatedly calling subroutine ssrkc4 for a single 
C classical fourth order RK step.  The truncation error is estimated 
C along the solution to adjust the integration step according to a 
C specified error tolerance.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   tf       final value of independent variable
C
C   u0       initial condition vector of length n
C
C   nsteps   maximum number of rkc4 steps
C
C   abserr   absolute error tolerance
C
C   relerr   relative error tolerance
C
C   u        ODE solution vector of length n after
C            nsteps steps
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size the arrays
      parameter(neq=250)
      dimension   e(neq)
      dimension    u0(n), u(n)
C
C Integration step
      h=(tf-t0)/2.0d0
C
C Minimum allowable step
      hmin=(tf-t0)/dfloat(nsteps)
C
C Start integration
      t=t0
C
C While independent variable is less than the final
C value, continue the integration
      do while(t.le.tf*0.999d0)
C
C       If the next step along the solution will go past 
C       the final value of the independent variable, set 
C       the step to the remaining distance to the final 
C       value
        if((t+h).gt.tf)then
          h=tf-t
        end if
C
C       Single rkc4 step
        call ssrkc4(n,t0,u0,h,t,u,e)
C
C       Flag for the end of the integration
        nfin1=1
C
C       Check if any of the ODEs have violated the error
C       critreria
        do i=1,n
          if(dabs(e(i)).gt.dabs(u(i))*relerr+abserr)then
C
C           Error violation, so integration is not complete.
C           Reduce integration step because of error violation
C           and repeat integration from the base point
            nfin1=0
            h=h/2.0d0
C
C           If the current step is less than the minimum
C           allowable step, set the step to the minimum
C           allowable value and continue integration from
C           new base point
            if(h.lt.hmin)then
              h=hmin
              nfin1=1
            end if
            go to 1
          end if
        end do
C
C     If there is no error violation, continue the integration
C     from new base point
 1    if(nfin1.eq.1)then
        do i=1,n
          u0(i)=u(i)
        end do
        t0=t
C
C       Test if integration step can be increased
        do i=1,n
          if(dabs(e(i)).gt.(dabs(u(i))*relerr+abserr)/16.0d0)then
C
C           Integration step cannot be increased
            go to 2
          end if
        end do
C
C       Increase integration step
        h=h*2.0d0
C
C     Continue for no error violation (nfin1=1)
 2    end if
C
C Continue do while
      end do
      return
C
C End of rkc4b
      end
      subroutine ssrkc4(n,t0,u0,h,t,u,e)
C
C Subroutine ssrkc4 computes an ODE solution by the classical fourth
C order RK method for one step along the solution (by calls to derv 
C to define the ODE derivative vector).  It also estimates the 
C truncation error of the solution, and applies this estimate as 
C a correction to the solution vector.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   u0       initial condition vector of length n
C
C   h        integration step
C
C   t        independent variable
C
C   u        ODE solution vector of length n after
C            one rkc4 step
C   
C   e        estimate of truncation error of the solu-
C            tion vector
C
C Double precision coding is used
      implicit double precision(a-h,k,o-z)
C
C Size the arrays
      parameter(neq=250)
      dimension ut0(neq), ut(neq), u4(neq),
     +           k1(neq), k2(neq), k3(neq), k4(neq)
      dimension    u0(n),    u(n),    e(n)
C
C Derivative vector at initial (base) point
      call derv(n,t0,u0,ut0)
C
C k1, advance of dependent variable vector and
C independent variable for calculation of k2
      do i=1,n
        k1(i)=h*ut0(i)
        u(i)=u0(i)+0.5d0*k1(i)
      end do
      t=t0+0.5d0*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k2, advance of dependent variable vector and
C independent variable for calculation of k3
      do i=1,n
        k2(i)=h*ut(i)
        u(i)=u0(i)+0.5d0*k2(i)
      end do
      t=t0+0.5d0*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k3, advance of dependent variable vector and
C independent variable for calculation of k4
      do i=1,n
        k3(i)=h*ut(i)
        u(i)=u0(i)+k3(i)
      end do
      t=t0+h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k4, stepping
      do i=1,n
        k4(i)=h*ut(i)
C
C       Second order step
        u(i)=u0(i)+k2(i)
C
C       Fourth order step
        u4(i)=u0(i)+(1.0d0/6.0d0)*(k1(i)+2.0d0*k2(i)+2.0d0*k3(i)+k4(i))
      end do
      do i=1,n
C
C       Truncation error estimate
        e(i)=u4(i)-u(i)
C
C       Fourth order solution vector (from 2,4 RK pair)
        u(i)=u(i)+e(i)
      end do
      return
C
C End of ssrkc4
      end
      subroutine rkf45a(n,t0,tf,u0,nsteps,u)
C
C Subroutine rkf45a computes an ODE solution by a fixed step
C rkf45 method for a series of points along the solution by 
C repeatedly calling subroutine ssrkf45 for a single rkf45 step.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   tf       final value of independent variable
C
C   u0       initial condition vector of length n
C
C   nsteps   number of rkf45 steps
C
C   u        ODE solution vector of length n after
C            nsteps steps
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size arrays
      parameter(neq=250)
      dimension   e(neq)
      dimension    u0(n), u(n)
C
C Integration step
      h=(tf-t0)/dfloat(nsteps)
C
C nsteps rkf45 steps
      do j=1,nsteps
C
C       Single rkf45 step
        call ssrkf45(n,t0,u0,h,t,u,e)
C
C       Reset base point values for the next rkf45 step
        do i=1,n
          u0(i)=u(i)
        end do
        t0=t
C
C     Next rkf45 step
      end do
C
C nsteps rkf45 steps completed
      return
C
C End of rkf45a
      end
      subroutine rkf45b(n,t0,tf,u0,nsteps,abserr,relerr,u)
C
C Subroutine rkf45b computes an ODE solution by a variable step
C RK Fehlberg method for a series of points along the solution 
C by repeatedly calling subroutine ssrkf45 for a single RK Fehlberg
C step.  The truncation error is estimated along the solution 
C to adjust the integration step according to a specified error 
C tolerance.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   tf       final value of independent variable
C
C   u0       initial condition vector of length n
C
C   nsteps   maximum number of rkf45 steps
C
C   abserr   absolute error tolerance
C
C   relerr   relative error tolerance
C
C   u        ODE solution vector of length n after
C            nsteps steps
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size the arrays
      parameter(neq=250)
      dimension   e(neq)
      dimension    u0(n), u(n)
C
C Initial integration step
      h=(tf-t0)/2.0d0
C
C Minimum allowable step
      hmin=(tf-t0)/dfloat(nsteps)
C
C Start integration
      t=t0
C
C While independent variable is less than the final
C value, continue the integration
      do while(t.le.tf*0.999d0)
C
C       If the next step along the solution will go past 
C       the final value of the independent variable, set 
C       the step to the remaining distance to the final 
C       value
        if((t+h).gt.tf)then
          h=tf-t
        end if
C
C       Single rkf45 step
        call ssrkf45(n,t0,u0,h,t,u,e)
C
C       Flag for the end of the integration
        nfin1=1
C
C       Check if any of the ODEs have violated the error
C       criteria
        do i=1,n
          if(dabs(e(i)).gt.dabs(u(i))*relerr+abserr)then
C
C           Error violation, so integration is not complete.
C           Reduce integration step because of error violation
C           and repeat integration from the base point
            nfin1=0
            h=h/2.0d0
C
C           If the current step is less than the minimum
C           allowable step, set the step to the minimum
C           allowable value and continue integration from
C           new base point
            if(h.lt.hmin)then
              h=hmin
              nfin1=1
            end if
            go to 1
          end if
        end do
C
C     If there is no error violation, continue the integration
C     from new base point
 1    if(nfin1.eq.1)then
         do i=1,n
            u0(i)=u(i)
         end do
         t0=t
C
C       Test if integration step can be increased
        do i=1,n
          if(dabs(e(i)).gt.(dabs(u(i))*relerr+abserr)/32.0d0)then
C
C           Integration step cannot be increased
            go to 2
          end if
        end do
C
C       Increase integration step
        h=h*2.0d0
C
C     Continue for no error violation (nfin1=1)
 2    end if
C
C Continue do while
      end do
      return
C
C End of rkf45b
      end
      subroutine ssrkf45(n,t0,u0,h,t,u,e)
C
C Subroutine ssrkf45 computes an ODE solution by the RK Fehlberg 45
C method for one step along the solution (by calls to derv to define 
C the ODE derivative vector).  It also estimates the truncation error 
C of the solution, and applies this estimate as a correction to the 
C solution vector.
C
C Argument list
C
C   n        number of first order ODEs
C
C   t0       initial value of independent variable
C
C   u0       initial condition vector of length n
C
C   h        integration step
C
C   t        independent variable
C
C   u        ODE solution vector of length n after
C            one rkf45 step
C   
C   e        estimate of truncation error of the solu-
C            tion vector
C
C Double precision coding is used
      implicit double precision(a-h,k,o-z)
C
C Size the arrays
      parameter(neq=250)
      dimension ut0(neq), ut(neq), u5(neq), 
     +           k1(neq), k2(neq), k3(neq), 
     +           k4(neq), k5(neq), k6(neq)
      dimension    u0(n),    u(n),    e(n)
C
C Derivative vector at initial (base) point
      call derv(n,t0,u0,ut0)
C
C k1, advance of dependent variable vector and
C independent variable for calculation of k2
      do i=1,n
        k1(i)=h*ut0(i)
        u(i)=u0(i)+0.25d0*k1(i)
      end do
      t=t0+0.25d0*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k2, advance of dependent variable vector and
C independent variable for calculation of k3
      do i=1,n
        k2(i)=h*ut(i)
        u(i)=u0(i)+(3.0d0/32.0d0)*k1(i)
     +            +(9.0d0/32.0d0)*k2(i)
      end do
      t=t0+(3.0d0/8.0d0)*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k3, advance of dependent variable vector and
C independent variable for calculation of k4
      do i=1,n
        k3(i)=h*ut(i)
        u(i)=u0(i)+(1932.0d0/2197.0d0)*k1(i)
     +            -(7200.0d0/2197.0d0)*k2(i)
     +            +(7296.0d0/2197.0d0)*k3(i)
      end do
      t=t0+(12.0d0/13.0d0)*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k4, advance of dependent variable vector and
C independent variable for calculation of k5
      do i=1,n
        k4(i)=h*ut(i)
        u(i)=u0(i)+( 439.0d0/ 216.0d0)*k1(i)
     +            -(   8.0d0         )*k2(i)
     +            +(3680.0d0/ 513.0d0)*k3(i)
     +            -( 845.0d0/4104.0d0)*k4(i)
      end do
      t=t0+h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k5, advance of dependent variable vector and
C independent variable for calculation of k6
      do i=1,n
        k5(i)=h*ut(i)
        u(i)=u0(i)-(   8.0d0/  27.0d0)*k1(i)
     +            +(   2.0d0         )*k2(i)
     +            -(3544.0d0/2565.0d0)*k3(i)
     +            +(1859.0d0/4104.0d0)*k4(i)
     +            -(  11.0d0/  40.0d0)*k5(i)
      end do
      t=t0+0.5d0*h
C
C Derivative vector at new u, t
      call derv(n,t,u,ut)
C
C k6, stepping
      do i=1,n
        k6(i)=h*ut(i)
C
C       Fourth order step
        u(i)=u0(i)+(  25.0d0 /  216.0d0)*k1(i)
     +            +(1408.0d0 / 2565.0d0)*k3(i)
     +            +(2197.0d0 / 4104.0d0)*k4(i)
     +            -(   1.0d0 /    5.0d0)*k5(i)
C
C       Fifth order step
        u5(i)=u0(i)+(   16.0d0/  135.0d0)*k1(i)
     +             +( 6656.0d0/12825.0d0)*k3(i)
     +             +(28561.0d0/56430.0d0)*k4(i)
     +             -(    9.0d0/   50.0d0)*k5(i)
     +             +(    2.0d0/   55.0d0)*k6(i)
      end do
      do i=1,n
C
C       Truncation error estimate
        e(i)=u5(i)-u(i)
C
C       Fifth order solution vector (from 4,5 RK pair)
        u(i)=u(i)+e(i)
      end do
      t=t0+h
      return
C
C End of ssrkf45
      end
