I have a dynamical system with one boundary saddle point, and one unstable interior point, and I would like to detect existence of cycles. This is of course a hard problem if the cycle is very small, and I can only hope for partial answers. My first attempt is a program with two inputs:the dynamical system, and the time T to iterate--- see below. I I use a very crude idea: starting paths near all the fixed points, and enlarging the time hoping to observe a cycle. If nothing else, this will detect separatices -- see the posts Plotting separatrices for nonlinear system , How to plot the stable and unstable manifolds of a hyperbolic fixed point of a nonlinear system of differential equations? The example below revealed a big cycle, but I have others from the literature where this crude program finds nothing (and contradicts the literature). So, my question was: if you needed to check quickly many numerics offered in the literature without accompanying notebooks (the typical case, unfortunately), what macro would you write?
T = 380; dyn = {16 - s/10 - (i s)/(100 (1 + i/1000)), -((21 i)/50) - (
13 i)/(2 (13/2 + i)) + (i s)/(100 (1 + i/1000))}; vars = {s, i};
equi = Solve[Thread[dyn == 0 && s >= 0 && i >= 0], vars] // N; le =
Length[equi];
Jacob = Grad[dyn, vars];
varst = Through[vars[t]];
diff = D[varst, t] \[Minus] (dyn /. Thread[vars -> varst]);
Print["Fixed points:"]
Do[p[j] = vars /. equi[[j]];
Jac[j] = Jacob /. Thread[vars -> p[j]];
Print[p[j], " with eigenvalues", Eigenvalues[Jac[j]]];
e[j] = Text["E[j]", Offset[{10, 11}, p[j]]];
st[j] = {PointSize[Large], Style[Point[p[j]], Green]};
b[j] = p[j] - {.1, 0}; a[j] = p[j] + {.1, 0};
eqb[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == b[j]]];
eqa[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == a[j]]];
eqb[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == b[j]]];
eqa[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == a[j]]];
sob[j] = NDSolve[eqb[j], vars, {t, 0, T}];
pb[j] =
ParametricPlot[{s[t], i[t]} /. sob[j], {t, 0, T},
PlotStyle -> Green] /.
Line[x_] :> {Arrowheads[Table[.03, {5}]], Arrow[x]};
soa[j] = NDSolve[eqa[j], vars, {t, 0, T}];
pa[j] =
ParametricPlot[{s[t], i[t]} /. soa[j], {t, 0, T},
PlotStyle -> Red] /.
Line[x_] :> {Arrowheads[Table[.05, {5}]], Arrow[x]}, {j, le}]
X = Table[p[j][[1]], {j, le}]; Y = Table[p[j][[2]], {j, le}]; xm =
3 Min[X]/4; xM = 99 Max[X]/98; ym = 0; yM = 4 Max[Y]/3;
epi = Flatten[Table[{e[j], st[j]}, {j, le}]];
sp = StreamPlot[{dyn}, {s, xm, xM}, {i, ym, yM}, StreamPoints -> 500,
Epilog -> epi, Frame -> True,
FrameLabel -> {"s", "i"},
LabelStyle -> Directive[Black, Medium]];
Show[sp, Table[pb[j], {j, le}], Table[pa[j], {j, le}]]
But, I accept Chris's implied answer that this problem is too difficult for one macro, and the quickest way is to use his package, which seems to provide even the period, thanks @Chris K