Chapter 3, Section 3.7

Here we define a general function, Stability, that (i) calculates the location of the specified equilibrium point (given by ptid) for the given mass ratio μ, (ii) calculates the eigenvalues of the displaced system and (iii) calculates the solution describing the motion around the equilibrium point for a displacement (X0,Y0).  The function uses most of the code from Location (see Fig.3.7.nb) to calculate the position of the equilibrium point.

Stability[μ_, ptid_Integer, {X0_, Y0_}] :=
Module[{α, r1, r2, μ1=1-μ, μ2=μ, lageq},
If[ptid==1,y=0;
  α=(μ/(3*(1-μ)))^(1/3);
  r20=α-(1/3)*α^2-(1/9)*α^3;
  lageq=μ1*(1-r2-1/(1-r2)^2)-μ2*(r2-1/r2^2);
  r2=r2/.FindRoot[lageq==0,{r2,r20}];
  r1=1-r2;
  x=μ1-r2];
If[ptid==2,y=0;
  α=(μ/(3*(1-μ)))^(1/3);
  r20=α+(1/3)*α^2-(1/9)*α^3;
  lageq=μ1*(1+r2-1/(1+r2)^2)+μ2*(r2-1/r2^2);
  r2=r2/.FindRoot[lageq==0,{r2,r20}];
  r1=1+r2;
  x=μ1+r2];
If[ptid==3,y=0;
  β=-(7/12)*(μ2/μ1)+(7/12)*(μ2/μ1)^2;
  r10=1+β;
  lageq=μ1*(r1-1/r1^2)+μ2*(1+r1-1/(1+r1)^2);
  r1=r1/.FindRoot[lageq==0,{r1,r10}];
  r2=1+r1; x=-r1-μ2];
If[ptid==4,x=1/2-μ2;
  y=N[Sqrt[3]/2];
  r1=1;
  r2=1];
If[ptid==5,x=1/2-μ2;
    y=N[-Sqrt[3]/2];
    r1=1;
    r2=1];
Print["Mass ratio = ", μ,"\n"];
Print["Lagrangian Point ", ptid,"\n"];
Print["{x0, y0} = ", {x,y},"\n"];
If[ptid == 1 || ptid == 2 || ptid ==3,
   abar = (1-μ)/r1^3 + μ/r2^3;
   ceqq = λ^4 + (2-abar)*λ^2 + (1+abar-2*abar^2);
   soln = Chop[λ /. NSolve[ceqq == 0, {λ}]]];
If[ptid == 4 || ptid == 5,
   ceqq = λ^4 + λ^2 + (27/4)*μ2*(1-μ2);
   soln = Chop[λ /. NSolve[ceqq == 0, {λ}]]];
abar = μ1/r1^3 + μ2/r2^3;
cbar = 3*(μ1*(x+μ2)/r1^5 + μ2*(x-μ1)/r2^5)*y;
dbar = 3*(μ1*(x+μ2)^2/r1^5 + μ2*(x-μ1)^2/r2^5);
uxx = 1 - abar + dbar;
uxy = cbar;
Print["Eigenvalues :"];
Print[ColumnForm[soln],"\n"];
Do[s[i] = soln[[i]], {i,1,4}];
Do[b[i] = (s[i]^2 - uxx)*a[i]/(2*s[i] + uxy), {i,1,4}];
eq1 = a[1] + a[2] + a[3] + a[4] == X0;
eq2 = b[1] + b[2] + b[3] + b[4] == Y0;
eq3 = s[1]*a[1] + s[2]*a[2] + s[3]*a[3] + s[4]*a[4] == 0;
eq4 = s[1]*b[1] + s[2]*b[2] + s[3]*b[3] + s[4]*b[4] == 0;
sol = Flatten[Chop[NSolve[
         {eq1,eq2,eq3,eq4},{a[1], a[2], a[3], a[4]}]]];
asol = {a[1], a[2], a[3], a[4]} /. sol;
bsol = {b[1], b[2], b[3], b[4]} /. sol;
Print["\nValues of α bar: "];
Print[ColumnForm[asol]];
Print["\nValues of β bar: "];
Print[ColumnForm[bsol]];
xx[t_] := 
Chop[ComplexExpand[
 ({a[1], a[2], a[3], a[4]} . 
    {Exp[s[1]*t], Exp[s[2]*t], Exp[s[3]*t], Exp[s[4]*t]}) /.
sol]];
yy[t_] := 
Chop[ComplexExpand[
 ({b[1], b[2], b[3], b[4]} . 
    {Exp[s[1]*t], Exp[s[2]*t], Exp[s[3]*t], Exp[s[4]*t]}) /.
sol]];
Print["\nX(t) = "];
Print[xx[t]];
Print["\nY(t) = "];
Print[yy[t],"\n"]
]

We consider two examples corresponding to the solutions given in Sect.3.7.1 and Sect.3.7.2 respectively.  In the first we examine the stability of the L1 point for a mass ratio μ=0.01 and small displacement of (0.00001,0.00001) from the equilibrium point.

Stability[0.01, 1, {0.00001,0.00001}]

In the second example we examine the stability of the L4 point for a mass ratio μ=0.01 and small displacement of (0.00001,0.00001) from the equilibrium point.

Stability[0.01, 4, {0.00001,0.00001}]