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

Question

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]]

Answer 1

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}}]
]

Answer 2

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]}]

Answer 3

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

Just adopting Chris answer to same sort of question,

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]