Problem finding separatrices for ODE system

Question

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

Answer 1

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.