(* ASTR 5540 Mathematical Methods Fall 1998.  Answer to Problem Set 7.

At the Mathematica prompt, type
<< ans7.math

*)

(* n = number of equations dY/dx = F(x,Y).
   Should be enough to use n = 2; then true for n > 2 by induction. *)
n = 2;

(*
 * Define n-dimensional vector versions Y of y and F of f.
 *)
Y[x_] := Table[y[i][x], {i,n}];
F[x_, y_] := Table[f[i][x, y], {i,n}];
f[i_][x__, {y__}] := f[i][x, y];	(* remove {} from arguments of f *)

(*
 * 4th order Runge-Kutta procedure.
 *)
Y0 := Y[x0];
F0 := F[x0, Y0];
F1 := F[x0 + dx/2, Y0 + dx/2 F0];
F2 := F[x0 + dx/2, Y0 + dx/2 F1];
F3 := F[x0 + dx, Y0 + dx F2];
YRK := Y0 + dx (F0/6 + F1/3 + F2/3 + F3/6);	(* next Y according to RK *)

(*
 * True solution.
 *)
Derivative[1][y[i_]][x_] := f[i][x, Y[x]];	(* dY/dx = F(x,y) *)
(* Next line tells Mathematica what it should already know,
   that the j'th derivative of Y is the j'th derivative of Y. *)
Derivative[j_ /; j >= 2][y[i_]][x_] := D[y[i][x], {x,j}];
Ytrue := Y[x0+dx];				(* next Y, true solution *)

(*
 * Compare RK and true solutions for Y to 4th order in dx;
   answer should be 0.
 *)
Series[ YRK - Ytrue , {dx,0,4} ]