f:=(t,y)->exp(t)*sin(y);RK4:=proc(ystart,tstart,tstop,deltat)
local t,y0,y1,k1,k2,k3,k4;
global results;
t:=tstart;
y0:=ystart;
results:=[[tstart,ystart]];
while t < tstop do
k1:=evalf(f(t,y0));
k2:=evalf(f(t+deltat/2,y0+k1*deltat/2));
k3:=evalf(f(t+deltat/2,y0+k2*deltat/2));
k4:=evalf(f(t+deltat,y0+k3*deltat));
y1:=evalf(y0+(k1+2*k2+2*k3+k4)*deltat/6);
t:=evalf(t+deltat);
results:=[op(results),op([[t,y1]])];
y0:=y1;
end do;
end proc;RK4(5,0,1,.05);