restart; eq:=D(y)=(s->f(s,y(s))); yp:=y(xn-h)+2*h*f(xn,y(xn)); ynp1:=y(xn)+h/2*(f(xn,y(xn))+f(xn+h,yp)); r:=y(xn+h)-ynp1; subs(h=0,r); T1:=diff(r,h); dr:=eval(subs(eq,T1)); subs(h=0,dr); ddr:=eval(subs(eq,diff(dr,h))); subs(h=0,ddr); d3r:=eval(subs(eq,diff(ddr,h))): subs(h=0,d3r);