#March 6, 2014; Starting Numerical solutions of ODEs

Help:=proc(): print(`Euler1(f,x,y,x0,y0,h),  `): 
print(`Euler(f,x,y,x0,y0,x1,h)`):
print(`impEuler1(f,x,y,x0,y0,h) `):
print(`impEuler(f,x,y,x0,y0,x1,h) `):
end:

#Euler1(f,x,y,x0,y0,h);
Euler1:=proc(f,x,y,x0,y0,h) 
y0+h*subs({x=x0,y=y0},f):
end:

#Euler(f,x,y,x0,y0,x1,h): Euler's method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of h
Euler:=proc(f,x,y,x0,y0, x1, h) local L,xc,yc:
L:=[y0]:

xc:=x0:
yc:=y0:

while xc<x1 do

yc:=Euler1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end:

#impEuler1(f,x,y,x0,y0,h);
impEuler1:=proc(f,x,y,x0,y0,h) local ystar:

ystar:=y0+h*subs({x=x0,y=y0},f):
y0+h*(subs({x=x0,y=y0},f)+subs({x=x0+h,y=ystar},f))/2:


end:

#impEuler(f,x,y,x0,y0,x1,h): Euler's method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of h
impEuler:=proc(f,x,y,x0,y0, x1, h) local L,xc,yc:
L:=[y0]:

xc:=x0:
yc:=y0:

while xc<x1 do

yc:=impEuler1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end: