#ProjectPrograms; Frank Wagner
Help:=proc(): print(` Thx(f,x,y,x0,y0,x1,h), RatMidPtEst(f,x,y,x0,y0,x1,m,h0) `): end:

#Thx(f,x,y,x0,y0,x1,h) inputs a function f in the variables x and y
#and four numbers x0, y0, x1, and h and
#OUTPUTS the expression T(h , 2*x1) for the initial condition ODE
#y’=f(x,y), y(x0)=y0.
Thx:=proc(f,x,y,x0,y0,x1,h) local l,a,i,n,x2:

x2:=2*x1:

l:=(x2-x0)/h:

if type(l,integer)<>true then
 RETURN(FAIL):
fi:

a[0]:=x0:

for i from 1 to l do
 a[i]:=a[i-1]+h:
od:

n[0]:=y0:
n[1]:=y0+h/2*subs({x=x0,y=y0},f):

for i from 2 to l do
 n[i]:=n[i-2]+h*subs({x=a[i-1],y=n[i-1]},f):
od:

evalf((1/2)*(n[l]+n[l-1]+h/2*subs({x=a[l],y=n[l]},f))):


end:


#RatMidPtEst(f,x,y,x0,y0,x1,m,h0) inputs a function f
#in the variables x and y, three numbers x0, y0, and x1,
#a positive integer m, and a mesh-size h0 and
#OUTPUTS the rational midpoint extrapolation estimate
#of the value y(x1) for the initial condition ODE
#y’=f(x,y), y(x0)=y0
RatMidPtEst:=proc(f,x,y,x0,y0,x1,m,h0) local T,h,i,k:

for i from 1 to m do
 T[-1,i]:=0:
od:

h[0]:=h0:
for i from 1 to m do
 h[i]:=h[0]/(i+1):
od:

for i from 0 to m do
 T[0,i]:=Thx(f,x,y,x0,y0,x1,h[i]):
od:

for k from 1 to m do
 for i from 0 to m-k do
  T[k,i]:=T[k-1,i+1]+(T[k-1,i+1]-T[k-1,i])/((h[i]/h[i+k])^2*
          (1-((T[k-1,i+1]-T[k-1,i])/(T[k-1,i+1]-T[k-2,i+1])))-1):
 od:
od:

T[m,0]:

end: