# Problem finding separatrices for ODE system

I have a system of two ODEs for which I would like to plot the separatrices. I tried the very neat method given here by Michael E2. Unfortunately, that code doesn't work for my system since the function sepICS[] returns only three starting points instead of four. Unfortunately I don't really understand what's going on in sepICS[] and thus I can't fix it.
The algorithm starts with the following code to find the saddle:

paramL = {sm -> 1.5*24, sw -> 1.*24, c -> 0.65, halfLwts -> 10, f -> 0.1*24, v1 -> 0.2*24, v2 -> 0.01*24, kfis -> 10, kfus -> 0.01, sf -> 0.5, sd -> 1.1};
sys = {(sw*wt)/((-mt*v2 + wt*(f - v2))/(c*v1) + 1) - (kfis/(kfis + kfus*(mt*sf + wt)))*wt*(Log[2]/halfLwts), (sm*mt)/((-mt*v2 + wt*(f - v2))/(c*v1) + 1) -  (kfis/(kfis + kfus*sf*(mt + wt)))*mt*sd*(Log[2]/halfLwts)} /. paramL;
vars = {wt, mt};
equilibria = Solve[sys == {0, 0}, vars, Reals];
saddles = Pick[equilibria, Sign[Det[D[sys, {vars}]]] /. equilibria, -1]


which gives: {{wt -> 997.648, mt -> 0}}.
A bit later the algorithm then uses this function:

sepICS[p0_, eps_] := With[{p1 = p0 + eps*Norm[p0]*{Cos[t],  Sin[t]}}, p1 /. NSolve[Det[{p1 - p0, sys /. Thread[vars -> p1]}] == 0 && 0<= t < 2*Pi]];


to generate 4 starting points around the saddle point (that are later used by NDSolve[]). But the problem is that for my system the function only generates 3 points

sepICS[{997.6476, 0}, 10^-7]
{{997.648, 0.0000814314}, {997.648, 1.22177*10^-20}, {997.648, -0.0000837899}}


or even only a single one

sepICS[{997.648, 0}, 10^-7]
{{997.648, 1.22177*10^-20}}


I'm not showing the rest of the code from the original post, since the problem lies in sepICS. So the question is: 1) Why is sepICS not resulting in four points for some systems and 2) What is the logic behind sepICS, why are exactly those starting points the right ones to trace back the separatrices ?

Many thanks

Sorry, I don't have time for a real answer right now, but here are some ideas that might help you or someone else come up with a robust solution. If you plot the function being solved inside sepICS you can see that it has four roots:

sepICS[p0_, eps_] := With[{p1 = p0 + eps*Norm[p0]*{Cos[t], Sin[t]}},
Print[Plot[
Det[{p1 - p0, sys /. Thread[vars -> p1]}], {t, 0, 2 π}]];
sol = NSolve[
Det[{p1 - p0, sys /. Thread[vars -> p1]}] == 0 && -1 <= t < 2 π];
Print[sol];
p1 /. sol];

sepICS[{997.6476, 0}, 10^-7]


Apparently NSolve is missing the root at or near t==0. I tried adjusting the constraints in NSolve but with no success. You might try another way to find all roots as discussed here.

Alternatively, and I'm not sure I understand @MichaelE2's code, could you get the initial conditions using Eigenvectors at the saddle point?

{λ, v} = Eigensystem[D[sys, {vars}] /. equilibria[[3]]]
(* {{-0.0173194, 0.00117657}, {{1., 0.}, {0.560376, 0.828239}}} *)


Here the first eigenvector would be the stable direction and the second the unstable direction relative to the equilibrium.

• (1) A simpler fix is to change the inequality in NSolve to 0 < t <= 2 Pi. I'm not sure either why 0 is rejected by NSolve. I suspect that there is rounding error in checking the interval, but I couldn't get it to go away even with WorkingPrecision -> 500. (2) I think, IIRC, the eigenvector method was too inaccurate in the old Q&A and produced a bad plot. Basically it represents an Euler-method step from the equilibrium along the separatrix, and my way is an Implicit-Euler step. I doubt that the Implicit Euler method is always superior for this type of problem. – Michael E2 Jun 17 at 14:52
• @MichaelE2 Thanks, I'll percolate on (2) later! – Chris K Jun 17 at 15:41
• Plotting the function inside sepICS is very helpful to understand the problem. And after changing the limits for t, everything is now working. Thanks everybody. – Axel Jun 18 at 9:09