/*
Runge-Kutta method for the differential equation
  y' = f(x,y)
where  f = f(x,y)  is a two-variable function on  x  and  y.

[Theorem]  The error of the Runge-Kutta method is  O(h^4).

[Proof]

Denoting the partial derivatives of  f  by
  f10 = f_{x},   f01 = f_{y},
  f20 = f_{xx},  f11 = f_{xy},  f02 = f_{yy},
  f30 = f_{xxx}, f21 = f_{xxy}, f12 = f_{xyy}, f03 = f_{yyy},
and the stride by  h,  compare the expansions of the true value 
and the approximate value of  y(x + h)  modulo  h^5.
*/

printf(false,"expansion of  f(x+s,y+t) :");
f_expand(s,t):=f+f10*s+f01*t+f20*s^2/2+f11*s*t+f02*t^2/2+f30*s^3/6+f21*s^2*t/2+f12*s*t^2/2+f03*t^3/6;

printf(false,"f_hat = f(x + h / 2, y + h * f / 2) :");
f_hat:remainder(expand(f_expand(h/2,h*f/2)),h^4);

printf(false,"f_tilde = f(x + h / 2, y + h * f_hat / 2) :");
f_tilde:remainder(expand(f_expand(h/2,h*remainder((f_hat/2),h^3))),h^4);

printf(false,"f_bar = f(x + h, y + h * f_tilde) :");
f_bar:remainder(expand(f_expand(h,h*remainder(f_tilde,h^3))),h^4);

printf(false,"The approximate value of  y(x + h) :");
y_approximate:expand(y+h*(f+2*f_hat+2*f_tilde+f_bar)/6);

printf(false,"y':");
y1:f;

printf(false,"y'':");
y2:f10+f01*f;

printf(false,"y''':");
y3:(f20+f10*f01)+(2*f11+f01^2)*f+f02*f^2;

printf(false,"y'''':");
y4:(f30+f20*f01+3*f10*f11+f10*f01^2)+(3*f21+5*f11*f01+3*f10*f02+f01^3)*f+(3*f12+4*f01*f02)*f^2+f03*f^3;

printf(false,"The true value of  y(x + h) :");
y_true:y+y1*h+y2*h^2/2+y3*h^3/6+y4*h^4/24;

printf(false,"The error value :");
expand(y_true-y_approximate);