#Execute read(`euler.txt`) from Maple and then Help(); to view the 
#procedures. f is an expression in x and y, x0 and y0 specify the
#initial condition, h is the step size, and xf is the final x-value.
#For example, try Euler(y^2,x,y,0,1,.5,1);

Help:=proc():
print(`Euler(f,x,y,x0,y0,h,xf),BackwardEuler(f,x,y,x0,y0,h,xf),RK4(f,x,y,x0,y0,h,xf) `):
end:

Euler:=proc(f,x,y,x0,y0,h,xf) local i,X,Y:
X:=[x0]:
Y:=[y0]:
for i from 1 while X[-1]<xf do
	X:=[op(X),X[-1]+h]:
	Y:=[op(Y),Y[-1]+h*subs({x=X[-2],y=Y[-1]},f)]:
od:
[X,Y]:
end:


BackwardEuler:=proc(f,x,y,x0,y0,h,xf) local i,X,Y,xn,yn:
X:=[x0]:
Y:=[y0]:
for i from 1 while X[-1]<xf do
	X:=[op(X),X[-1]+h]:
	Y:=[op(Y),solve(z=Y[-1]+h*subs({x=X[-1],y=z},f),z)]:
od:
[X,Y]:
end:

RK:=proc(f,x,y,x0,y0,h,xf) local i,X,Y,xn,yn,k1,k2,k3,k4:
X:=[x0]:
Y:=[y0]:
for i from 1 while X[-1]<xf do
	xn:=X[-1]:
	yn:=Y[-1]:
	k1:=subs({x=xn,y=yn},f):
	k2:=subs({x=xn+h/2,y=yn+h*k1/2},f):
	k3:=subs({x=xn+h/2,y=yn+h*k2/2},f):
	k4:=subs({x=xn+h/2,y=yn+h*k2/2},f):
	k4:=subs({x=xn+h,y=yn+h*k3},f):

	X:=[op(X),xn+h]:
	Y:=[op(Y),yn+h*(k1+2*k2+2*k3+k4)/6]:
od:
[X,Y]:
end: