(* traptrunc.wl -- Verify the trapezoid method is order 2 *)

eq=y'->Function[s,f[s,y[s]]]
r=y[t+h]-y[t]-h/2*(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. *)

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

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

Simplify[d3r/.h->0]
