ssrkf45:=proc(neqn,t0,u0,h,u,e)
#
# Procedure ssrkf45 computes an ODE solution by the RK Fehlberg 45
# 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
#
#   t        independent variable
#
#   u        ODE solution vector of length neqn after
#            one rkf45 step
#   
#   e        estimate of truncation error of the solu-
#            tion vector
#
# Type variables
  local ut0, ut, t, i, k1, k2, k3, k4, k5, k6, sum4, sum5:
#
# Declare arrays
  ut0:=array(1..neqn):   ut:=array(1..neqn):
   k1:=array(1..neqn):   k2:=array(1..neqn):
   k3:=array(1..neqn):   k4:=array(1..neqn):
   k5:=array(1..neqn):   k6:=array(1..neqn):
 sum4:=array(1..neqn): sum5:=array(1..neqn):
#
# Derivative vector at initial (base) point
  read "c:\\odelib\\maple\\pdenon\\derv.txt":
  derv(neqn,t0,u0,ut0):
#
# k1, advance of dependent variable vector and
# independent variable for calculation of k2
  for i from 1 to neqn do
    k1[i]:=h*ut0[i]:
     u[i]:=u0[i]+0.25*k1[i]:
  end do:
  t:=t0+0.25*h:
#
# Derivative vector at new u, t
  derv(neqn,t,u,ut):
#
# k2, advance of dependent variable vector and
# independent variable for calculation of k3
  for i from 1 to neqn do
    k2[i]:=h*ut[i]:
     u[i]:=u0[i]+(3.0/32.0)*k1[i]\
                +(9.0/32.0)*k2[i]:
  end do:
  t:=t0+(3.0/8.0)*h:
#
# Derivative vector at new u, t
  derv(neqn,t,u,ut):
#
# k3, advance of dependent variable vector and
# independent variable for calculation of k4
  for i from 1 to neqn do
    k3[i]:=h*ut[i]:
     u[i]:=u0[i]+(1932.0/2197.0)*k1[i]\
                -(7200.0/2197.0)*k2[i]\
                +(7296.0/2197.0)*k3[i]:
  end do:
  t:=t0+(12.00/13.0)*h:
#
# Derivative vector at new u, t
  derv(neqn,t,u,ut):
#
# k4, advance of dependent variable vector and
# independent variable for calculation of k5
  for i from 1 to neqn do
    k4[i]:=h*ut[i]:
     u[i]:=u0[i]+( 439.0/ 216.0)*k1[i]\
                -(   8.0       )*k2[i]\
                +(3680.0/ 513.0)*k3[i]\
                -( 845.0/4104.0)*k4[i]:
  end do:
  t:=t0+h:
#
# Derivative vector at new u, t
  derv(neqn,t,u,ut):
#
# k5, advance of dependent variable vector and
# independent variable for calculation of k6
  for i from 1 to neqn do
    k5[i]:=h*ut[i]:
     u[i]:=u0[i]-(   8.0/  27.0)*k1[i]\
                +(   2.0       )*k2[i]\
                -(3544.0/2565.0)*k3[i]\
                +(1859.0/4104.0)*k4[i]\
                -(  11.0/  40.0)*k5[i]:
  end do:
  t:=t0+0.5*h:
#
# Derivative vector at new u, t
  derv(neqn,t,u,ut):
#
# k6
  for i from 1 to neqn do
    k6[i]:=h*ut[i]:
  end do:
#
# Fourth order step
  for i from 1 to neqn do
    sum4[i]:=u0[i]+(    25.0/ 216.0)*k1[i]\
                  +(  1408.0/2565.0)*k3[i]\
                  +(  2197.0/4104.0)*k4[i]\
                  -(     1.0/   5.0)*k5[i]:
  end do:
#
# Fifth order step
  for i from 1 to neqn do
    sum5[i]:=u0[i]+(   16.0/  135.0)*k1[i]\
                  +( 6656.0/12825.0)*k3[i]\
                  +(28561.0/56430.0)*k4[i]\
                  -(    9.0/   50.0)*k5[i]\
                  +(    2.0/   55.0)*k6[i]:
  end do:
  t:=t0+h;
#
# Truncation error estimate
  for i from 1 to neqn do
    e[i]:=sum5[i]-sum4[i]:
  end do:
#
# Fifth order solution vector (from 4,5 RK pair);
# two ways to the same result are listed
  for i from 1 to neqn do
#   u[i]:=sum4[i]+e[i]:
    u[i]:=sum5[i]:
  end do:
#
# End of ssrkf45
  end: