(* Question 3 *)
eq=y'-> Function[t,f[t,y[t]]]
k1=f[tn,y[tn]]
k2=f[tn+h/3,y[tn]+h*k1/3]
k3=f[tn+2*h/3,y[tn]+2*h*k2/3]
ynp1=y[tn]+h*(k1+3*k3)/4
tau=y[tn+h]-ynp1

(* check the method is consistent *)
tau/.h->0

(* now take derivatives until the corresponding
  terms in the Taylor series cease to vanish *)
dj=tau;
For[j=1,j<6,j++,
dj=D[dj,h]/.eq;
djh0=Simplify[dj/.h->0];
Print["d",j,"h0=",djh0];
If[djh0=!=0,
Print["truncation error is O(h^",j,")"];
Print["order of method is O(h^",j-1,")"];
Break[]]]