sseuler:=proc(neqn,t0,u0,h,u,e)
#
# Procedure sseuler computes an ODE solution by the modified Euler 
# method for one step along the solution (by calls to derv to define 
# the ODE derivative vector).  It also estimates the truncation error 
# of the solution, and applies this estimate as a correction to the 
# solution vector.
#
# Argument list
#
#   neqn     number of first order ODEs
#
#   t0       initial value of independent variable
#
#   u0       initial condition vector of length neqn
#
#   h        integration step
#
#   u        ODE solution vector of length neqn after
#            one modified Euler step
#   
#   e        estimate of truncation error of the solu-
#            tion vector
#
# Type variables
  local ut0, ut, t, i:
#
# Declare arrays
  ut0:=array(1..neqn): ut:=array(1..neqn):
#
# Derivative vector at initial (base) point
  read "c:\\odelib\\maple\\pdenon\\derv.txt":
  derv(neqn,t0,u0,ut0):
#
# First order (Euler) step
  for i from 1 to neqn do
    u[i]:=u0[i]+ut0[i]*h:
  end do:
  t:=t0+h:
#
# Derivative vector at advance point
  derv(neqn,t,u,ut):
#
# Truncation error estimate
  for i from 1 to neqn do
    e[i]:=(ut[i]-ut0[i])*h/2.0:
  end do:
#
# Second order (modified Euler) solution vector
  for i from 1 to neqn do
    u[i]:=u[i]+e[i]:
  end do:
#
# End of sseuler
  end: