(* Stability of ODE's and methods *)

(* To run on command line:  math -run "<<ODE.stability.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`


(* An ODE is said to be stable over a certain interval if two solutions starting
 out close to each other remain close over the entire interval *)


(* Main case to explore: x' = a*x  *)  
(* To explore an unstable case, set a=1. Set x[0] = 1 or 1.1. For the stable
case set a = -1, and try x[0] = 1 and 1.1 *) 

a = -1;


solutionA = NDSolve[{x'[t] == a*x[t], x[0] == 1.0}, x, {t, 0, 10}]

plotA = Plot[Evaluate[x[t] /. solutionA ], {t, 0, 10}, PlotStyle->{Hue[1], Thickness[0.01]}, PlotLabel -> "x(0)= 1.0",  AxesLabel -> {t , x[t]}]

solutionB = NDSolve[{x'[t] == a*x[t], x[0] == 1.1}, x, {t, 0, 10}]

plotB = Plot[Evaluate[x[t] /. solutionB ], {t, 0, 10}, PlotStyle->{Hue[0.5], Thickness[0.01]}, PlotLabel -> "x(0) = 1.1",  AxesLabel -> {t, x[t]}]

(* Print the difference between the solutions at t = 10 *)

Print[]
Print["Difference between solutions at t=10     ", Evaluate[x[10] /. solutionB ] - Evaluate[x[10] /. solutionA ] ]
Print[]

Plot[Evaluate[x[t] /. solutionB ] - Evaluate[x[t] /. solutionA ], {t, 0, 10}, PlotLabel -> "Difference", AxesLabel -> {t, x[t]}]

Clear[a]