(***************************************************************************
Mathematica & C in modern Theoretical Physics
File:    SimpleNDSolve.m
Author:  J.-P. Kuska
Date:    /01/97
Version: Mathematica 2.2, 3.0
****************************************************************************)
BeginPackage["SimpleNDSolve`"]

EulerStep::usage="EulerStep[f_,y_,y0_,t_,tn_,h_] integrates over the
time interval tin [tn,tn+h] the systems of equations y'(t)==f,
with inital conditions y0 with the explicit Euler method."

RungeStep::usage="RungeStep[f_,y_,y0_,t_,tn_,h_] integrates over the
time interval tin [tn,tn+h] the systems of equations y'(t)==f,
with inital conditions y0 with the explicit Runge's method."

RungeKuttaStep::usage="RungeKuttaStep[f_,y_,y0_,t_,tn_,h_] integrates over the
time interval tin [tn,tn+h] the systems of equations y'(t)==f,
with inital conditions y0 with the explicit classical Runge-Kutta method."

SimpleNDSolve::usage="SimpleNDSolve[eqn,y,t,h,n] solves for
the differential equations eqn, depending on y[t] and t in
the interval t in [t0,t0+n*h] with constant step size h using
the Euler method. SimpleNDSolve[eqn,y,t,h,n,Method]
do the same as before but uses Method (RungeStep) instead of
the Euler method. The initial conditions must be given as
for Mathematicas NDSolve. SimpleNDSolve
returns a list {{t0,y[1][t0],...},{t1,y[1][t1],...},...}."

Begin["`Private`"]

(* get explicit differential equations *)
PrepareDGL[deqn:{___Equal},y_List,t_]:=
  Module[{dy},
         dy=D[#,t] & /@ y;
         dy //. Solve[deqn,dy]
    ] /; Length[deqn]==Length[y]

(* seperate initial conditions form the differential equations *)
SplitInitial[eqn:{___Equal},y_List,t_]:=
  Module[{ini=eqn,deqn,ydy},
         ydy=Join[y,D[#,t] & /@ y];
         While[ydy=!={},
           ini=Select[ini,FreeQ[#,First[ydy]] &];
           ydy=Rest[ydy];
          ];
         deqn=Complement[eqn,ini];
         {ini,deqn}
    ]

(* determine the starting time t0 *)
GetArgument[expr_,y_List,t_]:=Flatten[GetArgument[expr,#,t] & /@ y]
GetArgument[expr_List,y_,t_]:=Flatten[GetArgument[#,y,t] & /@ expr]
GetArgument[a_==b_,y_,t_]:={GetArgument[a,y,t],
                            GetArgument[b,y,t]}
GetArgument[expr_,y_,t_]:=
  Module[{tt,z},
         z=y /. t -> tt_;
         If[MatchQ[expr, z],
              First[expr],
              {}
           ]
        ]

(* Euler's method *)
EulerStep[f_,y_,y0_,t_,tn_,h_]:=
       y0+h* N[f /. Append[Thread[y->y0],t->tn] ]

(* Runge's method *)
RungeStep[f_,y_,y0_,t_,tn_,h_]:=
  Module[{k1,k2},
          k1=N[f /. Append[Thread[y->y0],t->tn]];
          k2=N[f /. Append[Thread[y->y0+h*k1],t->tn+h]];
          y0+0.5*h*(k1+k2)
        ]

RungeKuttaStep[f_, y_, y0_, t_, tn_, h_] := 
  Module[{k1, k2, k3, k4},
    k1 = N[f /. Append[Thread[y -> y0], t -> tn]];
    k2 = N[f /. Append[Thread[y -> y0 + h*k1/2], t -> tn + h/2]];
    k3 = N[f /. Append[Thread[y -> y0 + h*k2/2], t -> tn + h/2]];
    k4 = N[f /. Append[Thread[y -> y0 + h*k3], t -> tn + h]];
    y0 + h*((k1 + k4)/6 + (k2 + k3)/3)
   ]

SimpleNDSolve::badini="Insufficient number of initial conditions.";
SimpleNDSolve::badeqn="Insufficient number of differential equations.";
SimpleNDSolve::badt0="Initial conditions must be given at the same time.";
SimpleNDSolve::multisol="Multiple solutions for the right side of the
differential equations."

(* driver for the numerical solution of inital value problems *)

SimpleNDSolve[eqn:{___Equal},y_List,
              t_,h_Real,n_Integer,Method_:EulerStep]:=
  Module[{ini,deqn,f,args,t0,y0,res,moresol},
         {ini,deqn}=SplitInitial[eqn,y,t];
         If[Length[y]!=Length[ini],
            Message[SimpleNDSolve::badini];
            Return[$Failed]
           ];
         If[Length[y]!=Length[deqn],
            Message[SimpleNDSolve::badeqn];
            Return[$Failed]
           ];
         f=PrepareDGL[deqn,y,t];
         If[Length[f]===1, (* then *)
            f=Flatten[f],
            Message[SimpleNDSolve::multisol];
            Return[$Failed];
           ];
         args=GetArgument[ini,y,t];
         If[Max[args]=!=Min[args],
            Message[SimpleNDSolve::badt0];
            Return[$Failed];
           ];
         t0=First[args];
         y0=Flatten[(y /. t->t0)
              /. Solve[ini,y /. t->t0]] //N;
         res={Prepend[y0,t0]};
         For[i=0, i<n, i++,
             y0=Method[f,y,y0,t,t0+i*h,h];
             AppendTo[res,Prepend[y0,t0+(i+1)*h]];
            ];
         res
        ]
End[]
EndPackage[]
Null