(* thetatrunc.wl -- Find the truncation error as a function of theta *)

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

(* Taylor-series expand tau=r around h=0 *)

r/.h->0
dr=D[r,h]/.eq
T2=dr/.h->0
Simplify[T2]
d2r=D[dr,h]/.eq
T3=d2r/.h->0
Simplify[T3]
