(* rkgctrunc.wl -- RK method based on 2-point Gauss Quadrature *)

eq=y'->Function[s,f[s,y[s]]]
k1=f[t,y[t]]
k2=f[t+c2*h,y[t]+c2*h*k1]
k3=f[t+c2*h,y[t]+c2*h/2*(k1+k2)]
k4=f[t+c3*h,y[t]+c3*h*k1]
k5=f[t+c3*h,y[t]+c3*h/2*(k1+k4)]
c2=(-1/Sqrt[3]+1)/2
c3=(1/Sqrt[3]+1)/2
r=y[t+h]-y[t]-h/2*(k3+k5)

(* 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]
