(* f on [0.4, 0.75] defined by Appendix 2C *)
f[t_]:= 81600(30 t - 23)(5t - 2)(3 - 4t)/1127;
(* Verify that this agrees with Appendix 2B *) Expand[f[t]]

(* Calculate f[t+h] *) f[t+h]

(* Calculate Diff(f,[t,t+h] *) f[t+h] - f[t]

(* Simplify the above expression *) Expand[f[t+h]-f[t]]

(* Calculate DQ(f,[t,t+h] *) DQ = Expand[(f[t+h]-f[t])/h]

(* Extract the leading term *) fprime = Coefficient[DQ,h,0]

(* Find the error term *) epsilon = DQ - fprime