eq=y'-> Function[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
r/.h->0
T1=D[r,h]
dr=T1/.eq
dr/.h->0
ddr=D[dr,h]/.eq
ddr/.h->0
d3r=D[ddr,h]/.eq;
d3r/.h->0
