Reference : Jacobian from Pokorny' s paper

This is simply the Jacobian and its eigenvalues from Pokorny' s paper.  The output is currently hidden, simply remove the semicolon after the appropriate line to view.

a=.;
p=.;
t=.;
jac[a_,p_,t_]:={{0,1,0,0,0,0},{(1/(Abs[-1-p+a Sin[t]]))-1,0,0,0,0,0},
{0,0,0,1,0,0},{0,0,(1/(Abs[-1-p+a Sin[t]]))-1,0,0,0},{0,0,0,0,0,1},{0,0,0,0,-1,0}};
MatrixForm[jac[a,p,t]];
Eigenvalues[jac[a,p,t]];

Since jac[a,p,t].M[t] gives a block diagonal matrix (with the same entries for the first and second block), we can integrate one of these 2x2 blocks instead (of dealing with a 6x6 matrix).  We can also ignore the last block, which is solvable by hand and does not give eigenvalues greater than 1 (by modulus).

This is just double-checking my matrix algebra.  Note that in Mathematica, there are two kinds of matrix multiplication.  Multiplying matrices with a “.” gives the usual matrix multiplication.  Multiplying with a “*” or a space gives termwise multiplication (top left times top left, and so on)

a=.;
p=.;
fp={{0,1},{(1/(Abs[-1-p+a Sin[t]]))-1,0}};
mat={{moo,mot},{mto,mtt}};
fp . mat//MatrixForm;

Our code starts here

This code behaves as follows:
For given values of ‘a’ and ‘p’ (we step through these in steps of ‘ares’ and ‘pres’ respectively),
Use the NDSolve command (Numerical Differential Equation Solve), we find the solutons of the system 
M’[t]=jac[t]*M[t]
 M[0]=IdentityMatrix
By default, the output of NDSolve is an interpolating function (or in the case of a system, as we have, four interpolating functions).  Using the /. command, we evaluate these functions at T=2 π and append them to a 2x2 matrix appropriately.  We then find the eigenvalues of this matrix using the Eigenvalue command (this could be done more efficiently, as there is a convenient result for the eigenvalues of a 2x2 matrix.  However, this method is more intuitive and does not significantly add to the computation time (which can actually be found by putting the entire code inside a Timing[_____] command).
We then append to a table the values ‘a’, ‘p’ and the max of the eigenvalues (as we discussed, the system is unstable only if Abs[λ]>1 for one of the eigenvalues, so we only need be concerned with the eigenvalue of the highest modulus).

Note: when deciding if a given {a,p} is unstable, we must remember that we are doing this numerically, and therefore eigenvalues that should be exactly 1 (in modulus) may actually be very slightly less/greater (turns out, for this code, they ALWAYS turn out greater, with the exception of a=p=0, giving all eigenvalues 1), so we must give a little leeway to λ.
To quickly reproduce the plot from Pokorny’s paper, we make a refined list, called ‘unstable’, containing the all sets of {a,p} that are unstable.  We then use ListPlot to see these graphically. 

To go a bit deeper than Pokorny’s paper, we can make a refined list, called ‘unstableevs’, containing all sets of {a. p, Max[Abs[λ]]}.  We can then look at contour plots, density plots, 3D plots, etc (these are in the next section.

pmax=1;
pres=.01;
ares=.01;
list={};
For[p=0,p≤pmax,p+=pres,{
For[a=0,a<1+p,a+=ares,{
sol=NDSolve[
{
moo'[t]==mto[t],
mot'[t]==mtt[t],
mto'[t]==moo[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
mtt'[t]==mot[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
moo[0]==1,mot[0]==0,mto[0]==0,mtt[0]==1

},
{moo[t],mot[t],mto[t],mtt[t]},{t,0,2Pi}];
matM={{
(moo[t]/.sol/.t->2Pi)[[1]],
(mot[t]/.sol/.t->2Pi)[[1]]},
{(mto[t]/.sol/.t->2Pi)[[1]],
(mtt[t]/.sol/.t->2Pi)[[1]]}};
AppendTo[list,{p,a,Max[Abs[Eigenvalues[matM]]]}]
}]}]
stable={};
unstable={};
unstableevs={};
unstlist={};
res=.001;
For[i=1,i≤Length[list],i+=1,{
If[list[[i,3]]<=1+res,AppendTo[stable,{list[[i,1]],list[[i,2]]}],
AppendTo[unstable,{list[[i,1]],list[[i,2]]}]&&AppendTo[unstableevs,{list[[i,3]]}]&&AppendTo[unstlist,list[[i]]]
]
}]

ListPlot[unstable,PlotLabel->"A Quick Check of Unstable Regime"]
Show[%,Plot[1+x,{x,0,.5}]];(*this displays the same plot as before, with a line showing 1+p, the limit on a*)

Test : find the eigenvalues for a specific {a, p} to test method against the exact results obtained in Pokorny’s paper for a=0 and find the error  Mexact[t]-M[t]

This section of code also allows one to view the behavior of each element of M[t] graphically for given values of a and p.

(*atest=0;
ptest=1;
solu=NDSolve[
{
moo'[t]==mto[t],
mot'[t]==mtt[t],
mto'[t]==moo[t]( (1/(Abs[-1-ptest+ atest Sin[t]]))-1),
mtt'[t]==mot[t] ((1/(Abs[-1-ptest+atest Sin[t]]))-1),
moo[0]==1,mot[0]==0,mto[0]==0,mtt[0]==1

},{moo[t],mot[t],mto[t],mtt[t]},{t,0,2Pi}];
matM={{
(moo[t]/.solu/.t->2Pi)[[1]],
(mot[t]/.solu/.t->2Pi)[[1]]},
{(mto[t]/.solu/.t->2Pi)[[1]],
(mtt[t]/.solu/.t->2Pi)[[1]]}};
matM//MatrixForm;
Eigenvalues[matM];
Print["Determinant of matM = ",Det[matM]]
q=(ptest/(1+ptest));
mexact={{Cos[2Pi Sqrt[q]],(1/(Sqrt[q]))Sin[2Pi Sqrt[q]]},{-Sqrt[q]Sin[2 Pi Sqrt[q]],Cos[2Pi Sqrt[q]]}};(*this is the exact solution when a=0*)
Print["The error Mexact[2Pi]-M[2Pi] is given by  ",(matM-mexact)//MatrixForm]
frameval=True;
size=Large;
{{Plot[Evaluate[moo[t]/.solu],{t,0,2Pi},PlotRange->All,PlotLabel->{"moo[t] for {a,p}=" [{atest,ptest}]},Frame->frameval,ImageSize->size],
Plot[Evaluate[mot[t]/.solu],{t,0,2Pi},PlotRange->All,PlotLabel->"mot[t] for {a,p}=" [{atest,ptest}],Frame->frameval,ImageSize->size]},{
Plot[Evaluate[mto[t]/.solu],{t,0,2Pi},PlotRange->All,PlotLabel->"mto[t] for {a,p}=" [{atest,ptest}],Frame->frameval,ImageSize->size],
Plot[Evaluate[mtt[t]/.solu],{t,0,2Pi},PlotRange->All,PlotLabel->"mtt[t] for {a,p}=" [{atest,ptest}],Frame->frameval,ImageSize->size]}}//MatrixForm*)

Neat Plots

ListPointPlot3D[unstlist,AxesLabel->{"a","p","|λmax|"},ImageSize->Large,Filling->Bottom,PlotLabel->"Unstable regime, height at each point is |λ max| "]

ListDensityPlot[unstlist,PlotLabel->"Unstable Regime",AspectRatio->(1/2),PlotRange->{{0,p},{0,1+p}},ImageSize->Large]
ListContourPlot[unstlist,Contours->30,PlotLabel->"Unstable Regime",AspectRatio->(1/2),PlotRange->{{0,p},{0,1+p}},ImageSize->Large]

Applying our system

We want to find a range of masses for mb that will give us the unstable behavior for our system.
This DOESN’T WORK.  We did this backwards.  Instead, set a.

our parameters

ms=.022;
k=21.0;
lo=.223;

(*desired a*) a=.01;
res=.01;
pres=.01;
pmin=.;
For[p=0,p<1,p+=pres,{
sol=NDSolve[
{
moo'[t]==mto[t],
mot'[t]==mtt[t],
mto'[t]==moo[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
mtt'[t]==mot[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
moo[0]==1,mot[0]==0,mto[0]==0,mtt[0]==1

},
{moo[t],mot[t],mto[t],mtt[t]},{t,0,2Pi}];
matM={{
(moo[t]/.sol/.t->2Pi)[[1]],
(mot[t]/.sol/.t->2Pi)[[1]]},
{(mto[t]/.sol/.t->2Pi)[[1]],
(mtt[t]/.sol/.t->2Pi)[[1]]}};
If[Max[Abs[Eigenvalues[matM]]]≥1+res,pmin=p]
If[pmin>0,Break]
}]
pmin

For[p=pmin,p<1,p+=pres,{
sol=NDSolve[
{
moo'[t]==mto[t],
mot'[t]==mtt[t],
mto'[t]==moo[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
mtt'[t]==mot[t]( (1/(Abs[-1-p+a Sin[t]]))-1),
moo[0]==1,mot[0]==0,mto[0]==0,mtt[0]==1

},
{moo[t],mot[t],mto[t],mtt[t]},{t,0,2Pi}];
matM={{
(moo[t]/.sol/.t->2Pi)[[1]],
(mot[t]/.sol/.t->2Pi)[[1]]},
{(mto[t]/.sol/.t->2Pi)[[1]],
(mtt[t]/.sol/.t->2Pi)[[1]]}};
If[Max[Abs[Eigenvalues[matM]]]≤ 1+res,pmax=p]
If[pmax>0,Break]
}]
digits=3;
mbmin=Round[(10^digits) mb/.Solve[pmin==(((mb+(ms/2)))/(k lo))9.8,mb][[1,1]]]/(10^digits) 1.;
mbmax=Round[(10^digits)mb/.Solve[pmax==(((mb+(ms/2)))/(k lo))9.8,mb][[1,1]]]/(10^digits) 1.;
Print["The minimum mass mb that gives unstable behavior is ", mbmin," kg"]
Print["The maximum mass mb that gives unstable behavior is ", mbmax," kg"]
mbmax-mbmin