(* rk3trunc.wl -- Verify an RK3 is order 3 *)

eq=y'->Function[s,f[s,y[s]]]
k1=f[t,y[t]]
k2=f[t+2/3*h,y[t]+2/3*h*k1]
k3=f[t+2/3*h,y[t]+2/3*h*k2]
r=y[t+h]-y[t]-h*(1/4*k1+3/8*k2+3/8*k3)

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

r/.h->0
dr=D[r,h]/.eq;
dr/.h->0
d2r=D[dr,h]/.eq;
Simplify[d2r/.h->0]
d3r=D[d2r,h]/.eq;
Simplify[d3r/.h->0]
d4r=D[d3r,h]/.eq;
Simplify[d4r/.h->0]
