
(*********************************************************************

        Adapted from
        Roman E. Maeder: Programming in Mathematica,
        Second Edition, Addison-Wesley, 1991.

 *********************************************************************)


BeginPackage["RungeKutta`"]

RungeKutta::usage = "RungeKutta[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t1, dt}]
	numerically integrates the ei as functions of the yi with inital values ai.
	The integration proceeds in steps of dt from 0 to t1.
	RungeKutta[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t, t0, t1, dt}] integrates
	a time-dependent system from t0 to t1."

Begin["`Private`"]

RKStep[f_, y_, y0_, dt_] :=
	Module[{ k1, k2, k3, k4 },
		k1 = dt N[ f /. Thread[y -> y0] ];
		k2 = dt N[ f /. Thread[y -> y0 + k1/2] ];
		k3 = dt N[ f /. Thread[y -> y0 + k2/2] ];
		k4 = dt N[ f /. Thread[y -> y0 + k3] ];
		y0 + (k1 + 2 k2 + 2 k3 + k4)/6
	]

RungeKutta[f_List, y_List, y0_List, {t1_, dt_}] :=
	NestList[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ] /;
		Length[f] == Length[y] == Length[y0]

RungeKutta[f_List, y_List, y0_List, {t_, t0_, t1_, dt_}] :=
	Module[{res},
		res = RungeKutta[ Append[f, 1], Append[y, t], Append[y0, t0], {t1 - t0, dt} ];
		Drop[#, -1]& /@ res
	]  /;  Length[f] == Length[y] == Length[y0]

End[]

Protect[RungeKutta]

EndPackage[]
