#rungekutta(f,t0,y0,h,n) computes approximate values to dy/dt=f, y(t0)=y0
#using the runge kutta method with step size h, for n iterations
#
#for example, to do #2 pn p. 461 where dy/dt=5t-3sqrt(y) y(0)=2, with h=0.1,
#try rungekutta(5*t-3*sqrt(y),0,2,0.1,4);
rungekutta:=proc(f,t0,y0,h,n) local i, k1,k2,k3,k4,tprev,yprev,tcurr,ycurr:
tcurr:=t0:
ycurr:=y0:

for i from 1 to n do
tprev:=tcurr:
yprev:=ycurr:

k1:=evalf(subs({t=tprev,y=yprev},f)):
k2:=evalf(subs({t=tprev+h/2,y=yprev+h/2*k1},f)):
k3:=evalf(subs({t=tprev+h/2,y=yprev+h/2*k2},f)):
k4:=evalf(subs({t=tprev+h,y=yprev+h*k3},f)):

#compute new value of y by equations (2) and (3) on p. 458
ycurr:=evalf(yprev+h*(k1+2*k2+2*k3+k4)/6):
tcurr:=evalf(tprev+h):

print("when t =", tcurr, " y is approximately ", ycurr):

od:

end: