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