rkf45b:=proc(neqn,t0,tf,u0,nsteps,abserr,relerr,u)
#
# Procedure rkf45b computes an ODE solution by the variable
# step RK Fehlberg 45 method for a series of points along
# the solution by repeatedly calling procedure ssrkf45 for
# a single RK Fehlberg 45 step. The truncation error is
# estimated along the solution to adjust the integration
# step according to a specified error tolerance.
#
# Argument list
#
# neqn number of first order ODEs
#
# t0 initial value of independent variable
#
# tf final value of independent variable
#
# u0 initial condition vector of length neqn
#
# nsteps maximum number number of rkf45 steps
#
# abserr absolute error tolerance
#
# relerr relative error tolerance
#
# u ODE solution vector of length neqn after
# nsteps steps
#
# Type variables
local e, h, i, j, tm, t, hmin, nfin1:
#
# Size arrays
e:=array(1..neqn):
#
# Initial integration step
h:=(tf-t0)/2.0:
#
# Minimum allowable step
hmin:=(tf-t0)/nsteps:
#
# Start integration
tm:=t0:
#
# While independent variable is less than the final
# value, continue the integration
while(tm <= tf*0.999) do
#
# If the next step along the solution will go past
# the final value of the independent variable, set
# the step to the remaining distance to the final
# value
if(tm+h > tf) then h:=tf-tm: end if:
#
# Single rkf45 step
read "c:\\odelib\\maple\\ode1x1\\ssrkf45.txt":
ssrkf45(neqn,tm,u0,h,u,e):
t:=tm+h:
#
# Flag for the end of the integration
nfin1:=1:
#
# Check if any of the ODEs have violated the error
# criteria
for i from 1 to neqn do
if(abs(e[i]) > abs(u[i])*relerr+abserr) then
#
# Error violation, so integration is not complete.
# Reduce integration step because of error violation
# and repeat integration from the base point
h:=h/2.0: nfin1:=0: break:
end if:
end do:
#
# If the current step is less than the minimum
# allowable step, set the step to the minimum
# allowable value and continue integration from
# new base point
if(h < hmin) then h:=hmin: nfin1:=1: end if:
#
# If there is no error violation, continue the integration
# from new base point
if(nfin1 = 1) then
for i from 1 to neqn do
u0[i]:=u[i]:
end do:
tm:=t:
#
# Test if integration step can be increased
for i from 1 to neqn do
if(abs(e[i]) > (abs(u[i])*relerr+abserr)/32.0) then
#
# Integration step cannot be increased
nfin1:=0: break:
end if:
end do:
#
# Increase integration step
if(nfin1 = 1) then h:=h*2.0: end if:
#
# Continue for no error violation (nfin1=1)
end if:
#
# Continue while
end do:
#
# End of rkf45b
end: