(* Stability of ODE's and methods *)

(* To run on command line:  math -run "<<ODE.stability1.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. Otherwise, the ODE is unstable.  *)



(* An unstable ODE to explore x' = x^2  Seek solution in the range 0 <= t <
0.999999 *)  

solutionA = NDSolve[{x'[t] == (x[t]*x[t]), x[0] == 1.0}, x, {t,0, 0.9999}]; 



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

solutionB = NDSolve[{x'[t] == (x[t]*x[t]), x[0] == 1.01}, x, {t,0, 0.9999}]; 

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

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

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

Plot[Evaluate[x[t] /. solutionB ] - Evaluate[x[t] /.
solutionA ], {t, 0, 0.99}, PlotRange -> {{0, 1.1},{0,
15000}}, PlotLabel -> "Difference", AxesLabel -> {t, x[t]}, PlotStyle->{Hue[0.1], Thickness[0.01]} ]

Clear[a]