Stiffness:  What does it look like!

We have a variety of test problems of the form
	y'=f[y]
Some are stiff and some are not.  How would we tell which is which.  Well we should look at the linearized equation
	y'=f[y0]+Df[y0].(y-y0)
whose behavoir is governed by the eigenvalues of the matrix M=Df[y0].  Clearly this varies from point to point along the
solution curve.

We want to see what happens along some of our test problem solution curves.

p60 Test problem

TMax=50;
fp60[{t_,y1_,y2_}]:={1,
-16 y1 + 12 y2 + 16 Cos[t] - 13 Sin[t],
12 y1 - 9 y2 - 11 Cos[t] + 9 Sin[t]};
M[{t_,y1_,y2_}]= D[fp60[{t,y1,y2}],{{t,y1,y2}}];
ySol= y/.NDSolve[{
y'[t]==fp60[y[t]],
y[0]=={0,1,0}
},
y,
{t, 0, TMax}][[1]];
Plot[ySol[t][[2;;3]],{t,0,TMax}]

This one is kinda silly because the derivative is independent of y1 and y2 and thewigen values are purely real.

StiffTest[t_]:=Module[
{zs},
zs=Eigenvalues[M[ySol[t]]]
]
Plot[ StiffTest[t],{t, 0,TMax}]

Gravity p29 test problem

Clear[GravityP29]
GravityP29[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=100/8145, d1,d2},
d1=((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(-3/2));
d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(-3/2));
{y4,y5,y6,
2y5 + y1 - μ(y1+μ-1)d1 -(1-μ)(y1+μ)d2,
-2y4 + y2 - μ y2 d1 -(1-μ)y2 d2,
- μ y3 d1 -(1-μ)y3 d2
}
]
IC1 = SetPrecision[{0.994,0,0,0,-2.001585106370825224,0},64];
IC=IC1;
TMax=17.06521656015796; 
ySol=y/.NDSolve[{
y'[t]==GravityP29[y[t]],
y[0]==IC
},
y,
{t,0, TMax},
PrecisionGoal->16,
AccuracyGoal->16,
InterpolationOrder->16,
WorkingPrecision->64][[1]];
TStart=TMax/3;
H=TMax/3;
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax},
Epilog->{PointSize[0.02],Pink, 
Point[ySol[TStart][[{1,2}]]],
Point[ySol[TStart + H][[{1,2}]]]
}
]

Clear[G,H]
MGravityP29[{y1_,y2_,y3_,y4_,y5_,y6_}]= D[GravityP29[{y1,y2,y3,y4,y5,y6}],{{y1,y2,y3,y4,y5,y6}}];
G[t_]:=Module[
{zs},
zs=Sort[Eigenvalues[ MGravityP29[ySol[t]]],Re[#1]≤Re[#2]&];
zs/.{
z_Complex:>{Re[z],Im[z]},
z_Real:>{z,0}
}
]
ParametricPlot[
G[t],
{t,0.1, 0.9TMax},
AspectRatio->1]

Plot[
{
G[t][[1,1]],
G[t][[1,2]],
G[t][[-1,1]],
G[t][[-1,2]]
},
{t,0.2, 0.9TMax},
AspectRatio->1,
PlotRange->All]