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

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

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

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]