      subroutine euler2a(n,t0,tf,u0,nsteps,u)
C
C Subroutine euler2a computes an ODE solution by a fixed step
C modified Euler method for a series of points along the solution 
C by repeatedly calling subroutine sseuler for a single modified
C Euler 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 modified Euler 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
      dimension u0(n), u(n), e(n)
C
C Integration step
      h=(tf-t0)/dfloat(nsteps)
C
C nsteps modified Euler steps
      do j=1,nsteps
C
C       Modified Euler step
        call sseuler(n,t0,u0,h,t,u,e)
C
C       Reset base point values for the next modified 
C       Euler step
        do i=1,n
          u0(i)=u(i)
        end do
        t0=t
C
C     Next modified Euler step
      end do
C
C nsteps modified Euler steps completed
      return
C
C End of euler2a
      end
      subroutine euler2b(n,t0,tf,u0,nsteps,abserr,relerr,u)
C
C Subroutine euler2b computes an ODE solution by a variable step
C modified Euler method for a series of points along the solution 
C by repeatedly calling subroutine sseuler for a single modified 
C Euler 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 modified Euler 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)/8.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 modified Euler step
        call sseuler(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)/4.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 euler2b
      end
      subroutine sseuler(n,t0,u0,h,t,u,e)
C
C Subroutine sseuler computes an ODE solution by the modified Euler 
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 modified Euler 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,o-z)
C
C     Size the arrays
      parameter(neq=250)
      dimension  ut(neq),ut0(neq)
      dimension    u0(n),    u(n),     e(n)
C
C Derivative vector at initial (base) point
      call derv(n,t0,u0,ut0)
C
C First order (Euler) step
      do i=1,n
         u(i)=u0(i)+ut0(i)*h
      end do
      t=t0+h
C
C Derivative vector at advance point
      call derv(n,t,u,ut)
C
C Second order step
      do i=1,n
C
C        Truncation error estimate
         e(i)=(ut(i)-ut0(i))*h/2.0d0
C
C        Second order (modified Euler) solution vector 
         u(i)=u(i)+e(i)
      end do
      return
C
C End of sseuler
      end
