(* rk2trunc.wl -- Verify an RK2 method is order 2 *)

eq=y'->Function[s,f[s,y[s]]]
r=y[t+h]-y[t]-h/2*(f[t,y[t]]+f[t+h,y[t]+h*f[t,y[t]]])

(* Taylor-series expand the truncation error tau=r about h=0
   until we find the first non-zero term. *)

f[t+h,y[t]+h*f[t,y[t]]])
r/.h->0
dr=D[r,h]/.eq
f[t,y[t]]+y[t]])
dr/.h->0
d2r=D[dr,h]/.eq
d2r/.h->0
d3r=D[d2r,h]/.eq;
d3r/.h->0

(* Double check r'''(0)<>0 with Simplify *)

Simplify[d3r/.h->0]
