(* ODE *) 
(* To run on command line:  math -run "<<ODEUpdated.math" *)
(* The above line assumes that math kernel is aliased to "math" *)
(* For example, on my machine alias math='/Applications/Mathematica.app/Contents/MacOS/MathKernel' *) 
(* or simply use << ODE.math once math kernel is already running *) 

<<JavaGraphics`

Remove[x]
(* First Example *)
(* Find the solution of a first order "key" ODE *)
(* x'= kx   *)

k = -1;
x0 = 1.0;  (* Starting Point *)
max = 10.0;  (* max range *) 


(* Mathematica's Analytical solution *)  
analytsol = DSolve[{x'[t] == k x[t], x[0] == x0 }, x[t], t]

gr3 = Plot[Evaluate[x[t]  /. analytsol], {t, 0, 10}, PlotRange -> All,
PlotStyle -> {Green}, DisplayFunction ->Identity];

Show[gr3, DisplayFunction -> $DisplayFunction, PlotLabel->"Exact solution" ]

(* Mathematica's Numerical solution. Here we will use Euler's method *)
(* Step size h is controlled by StartingStepSize, do not change other 
parameters *)                           

numsol = NDSolve[{x'[t] == k x[t], x[0] == x0 }, x, {t, 0, max}, Method ->
{"FixedStep", Method -> "ExplicitEuler"}, StartingStepSize -> 0.5, MaxStepFraction -> 1 ]

 
gr4 = Plot[Evaluate[x[t] /. numsol], {t, 0, max}, PlotStyle -> {Red}, PlotRange -> All, DisplayFunction -> Identity];

Show[gr4, DisplayFunction -> $DisplayFunction, PlotLabel->"Numerical solution" ]

Show[gr3, gr4, DisplayFunction -> $DisplayFunction, PlotLabel->"Exact
(green) and Numerical (Red) solutions together" ]
 

(* Now hand code the simplest Euler's method *) 
h = 0.1; (* time discretization step *)
s= max/h;    (* number of steps: note that it grows as h gets smaller *)  
yh[0] = 1;(* initial value *)

Do[yh[n + 1] = yh[n] + h*(k yh[n]), {n, 0, s-1}]

Eulersol = Table[{N[h i], yh[i]}, {i, 0, s}] 

gr5 = ListPlot[Eulersol, PlotStyle -> {Red, PointSize[0.02]},
PlotLabel->"Solution by hand-coded Euler's method"]

Show[gr3, gr5, DisplayFunction -> $DisplayFunction, PlotLabel->"Exact
(green) vs.  hand coded Euler's method (red)", AxesLabel -> {Style["t", Bold, 16] , Style["y(t)", Bold, 16]}]

(* What happens if k > 0 ? *) 

(* Second Example *)
(* Find the solution of a second order "key" ODE *)
(* x’’=kx *)
(* For k< 0 it descibes a harmonic oscillator with frequency w=sqrt(|k|), 
and period  T= 2Pi/sqrt(|k|). x(t) = cos(wt) *)  


Remove[x]
k = -4;
max2 = 40; (* range *) 

analytsol2 = DSolve[{x''[t] == k x[t], x[0] == 1, x'[0] == 0}, x[t], t]

numsol2 = NDSolve[{x''[t] == k x[t], x[0] == 1, x'[0] == 0}, {x}, {t, max2}]

gr6=Plot[Evaluate[x[t]/.analytsol2],{t,0,max2},PlotRange->All,PlotStyle->{Green},DisplayFunction->Identity, PlotLabel->"exact solution"];

Show[gr6, DisplayFunction->$DisplayFunction]

(* 
gr7=Plot[Evaluate[x[t]/.numsol2],{t,0,max2},PlotRange->All,DisplayFunction->Identity, PlotLabel->"Mathematica's numerical solution" ];

Show[gr6,gr7,DisplayFunction->$DisplayFunction]
*) 

(* Hand coded Euler's *) 

h=0.08; 
s= max2/h; (* number of steps need to reach the same range max2 *) 
yh[0]=0; xh[0]=1; (* initial values *)

Do[ xh[n+1]= xh[n] + h*(yh[n]);yh[n+1]=yh[n]+h*(k xh[n]),{n,0,s-1}]

Eulersol2=Table[{N[h i],xh[i]},{i,0,s}]

gr8=ListPlot[Eulersol2,PlotStyle->{Red}, PlotLabel->"hand coded Euler's method"]


Show[gr6, gr8, DisplayFunction->$DisplayFunction, PlotLabel->"Exact (green)
vs. hand coded Euler's (red) ", AxesLabel -> {Style["t", Bold, 16] , Style["y(t)", Bold, 16]} ]

(*  save tables to file if needed
Export["numSol.data", Eulersol, "Table"]
Export["exactSol.data", 
 Table[{N[t], Evaluate[x[t] /. analytsol]}, {t, 0, 10}], "Table"]
*)