(* Stability of ODE's and methods *) 

(* To run on command line:  math -run "<<ODE.stability2.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 ru
nning *)

<<JavaGraphics`


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

(* An ODE can be stable in one region, and unstable in another *)

(* For example, the one below is unstable if t < 1, and stable if t > 1 *)

(* Even though a particular ODE {x' == f(t, x)} may be stable, a method applied to it may be unstable.   *)

solutionA = NDSolve[{x'[t] == x[t]*x[t] -t, x[-1.99] == 1.0}, x, {t, -2,
0}]

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


solutionB = NDSolve[{x'[t] == x[t]*x[t] -t, x[-1.99] == 1.00000001}, x, {t,
-2, 0}]

plotB = Plot[Evaluate[x[t] /. solutionB ], {t, -2, 0}, PlotStyle->{Hue[0.5],
Thickness[0.01]}, PlotLabel -> "x(-1.99) = 1.00000001",  AxesLabel -> {t, x[t]}, PlotRange -> All]

Show[plotA, plotB, PlotLabel -> "both solutions", AxesLabel -> {t, x[t]},
PlotRange -> All ] 


Plot[Evaluate[x[t] /. solutionB ] - Evaluate[x[t] /. solutionA ], {t, -2, 0}, PlotLabel -> "Difference between them",  AxesLabel -> {t, x[t]}, PlotRange -> All]