# How to plot the stable and unstable manifolds of a hyperbolic fixed point of a nonlinear system of differential equations?

Suppose we have the following simplified system of two ordinary differential equations:

$$\dot{x}(t)=x(t)^2+2y(t)\\ \dot{y}(t)=3x(t)$$

The system has a hyperbolic fixed point the origin. Hence there exits a stable and an unstable manifold going through the origin.

I would like to plot these exact sets. I am able to plot basins of attraction for regions in phase space but I fail to do so once the attracting region is a one-dimensional line (as is the case here with the stable manifold).

Edit: With MMM's help I was able to roughly plot the unstable manifold. However, what is more challenging is to plot the stable one.

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6}];
Show[splot,
ParametricPlot[
Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x,
y}, {t, 0, 3}]]], {t, 0, 3}, PlotStyle -> Red],
ParametricPlot[
Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x,
y}, {t, 0, 2.4}]]], {t, 0, 2.4}, PlotStyle -> Red]]


• Share your try in Mathematica.
– zhk
Commented Jan 13, 2017 at 10:15

The trick to getting the unstable manifold is to run backwards in time. Here I added two more solutions starting near the equilibrium but with negative t:

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6}];
Show[splot,
ParametricPlot[Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x, y}, {t, 0, 4}]]],
{t, 0, 4}, PlotStyle -> Red, PlotRange -> {{-6, 6}, {-6, 6}}],
ParametricPlot[Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x, y}, {t, 0, 2.6}]]],
{t, 0, 2.6}, PlotStyle -> Red, PlotRange -> {{-6, 6}, {-6, 6}}],
ParametricPlot[Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x, y}, {t, 0, -2.6}]]],
{t, 0, -2.6}, PlotStyle -> Green, PlotRange -> {{-6, 6}, {-6, 6}}],
ParametricPlot[Evaluate[First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x, y}, {t, 0, -4}]]],
{t, 0, -4}, PlotStyle -> Green, PlotRange -> {{-6, 6}, {-6, 6}}]
]


• Thank you, this is exactly what I wanted. Commented Jan 13, 2017 at 13:26

Here's one way to get StreamPlot to show the desired manifolds and fill in the rest of the plot around them.

We can use the eigenvectors of the Jacobian at the equilibrium to approximate a point on each manifold. From such a point, another point further away from the equilibrium along the manifold may be constructed with NDSolve. Using these for StreamPoints allows StreamPlot to fill in the rest of the phase portrait.

(* some useful quantities/expression for playing *)
sys = {x'[t] == x[t]^2 + 2*y[t], y'[t] == 3*x[t]};
vars = {x, y};
dvars = D[Through[vars[t]], t];                   (* {x'[t], y'[t]} *)
vel = dvars /. First@Solve[sys, dvars];           (* expressions for {x'[t], y'[t]} *)
equil = First@Solve[sys /. Thread[dvars -> 0]];   (* equilibrium *)
phasefield = vel /. var_[t] :> var;               (* strip [t] from expressions *)
linearized = D[vel, {Through[vars[t]]}] /. equil; (* its Eigensystem[] describes equil. *)

(* compute stream points for the separatrix manifolds *)
sp = Function[{λ,      (* eigenvalue *)
v,      (* eigenvector *)
side},  (* ±1 = which side of the equilibrium: ± v *)
Module[{tfinal},
{NDSolveValue[{sys, Through[vars[0]] == 10^-4 v*side,
WhenEvent[Norm[#] == 1, tfinal = t; "StopIntegration"] &@
Through[vars[t]]},
Through[vars[tfinal]], {t, 0, 100 Sign@λ}],
If[λ < 0, Green, Red]}
]
] @@@ Append @@@ Tuples[{Thread@Eigensystem@linearized, {-1, 1}}]
(*
{{{ 0.599966, -0.800025}, RGBColor[0, 1, 0]},
{{-0.664941,  0.746896}, RGBColor[0, 1, 0]},
{{-0.599966, -0.800025}, RGBColor[1, 0, 0]},
{{ 0.664941,  0.746896}, RGBColor[1, 0, 0]}}
*)

splot = StreamPlot[phasefield, {x, -6, 6}, {y, -6, 6},
StreamPoints -> {Append[sp, Automatic]}]


• Nice! I think there may be a typo or two in the code, though...it's not working for me. Commented Jan 15, 2017 at 6:01
• @SimonRochester Thanks for letting me know. I had tried to simplify a couple of things and forgot to change all references to old variables. Commented Jan 15, 2017 at 13:29

Your question is not clear. But here is a starting for you.

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6},
StreamColorFunction -> "Rainbow"]
Manipulate[
Show[splot,
ParametricPlot[
Evaluate[
First[{x[t], y[t]} /.
NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == point]}, {x, y}, {t,
0, T}]]], {t, 0, T}, PlotStyle -> Red]], {{T, 1}, 0.1,
1}, {{point, {0.01, 0.0}}, Locator}, SaveDefinitions -> True]


• more specifically: I want to plot the two (unique) trajectories that go through the saddle at the origin. One of them is the stable manifold, the other one is unstable. By studying the vector field we can get a good idea of how they look like, but i want to plot them. Commented Jan 13, 2017 at 11:27
• @ParetoWilli I think you want something exactly like this. mathematica.stackexchange.com/questions/80284/…
– zhk
Commented Jan 13, 2017 at 12:43