StreamPoints for all of the domain (how to color basins of attraction)

Horribly slow, but you probably get the idea. If this system can be solved analytically, then one could implement a far better approach, but I haven't looked into this.

What you are after is to find out which point in your plane converges to which sink. You can do this by solving the differential equation where your direction vector is the speed. You provide a starting point and start integrating. If you have done this, you can look what point is reached at the end of the integration and use the specific color for the sink.

A function that does exactly this is the following. It additionally stores the path starting from the initial point to the end.

line[{x0_, y0_}] :=
 Module[{x, y, t, res, col},
  res = NDSolveValue[{x'[t] == y[t], 
     y'[t] == -.15 y[t] + x[t] - x[t]^3, x[0] == x0, y[0] == y0}, {x, 
     y}, {t, 0, 100}];
  col = If[Norm[{-1, 0} - Through[res[100]]] < .01, Green, Red];
  {col, Line[Table[Through[res[t]], {t, 0, 100, .1}]]}
  ]

Now you can sample your region and for each starting point you create a colored line (this may take a while)

lines = ParallelTable[line[{x, y}], {x, -2, 2, .1}, {y, -2, 2, .1}];

After that, you can combine all lines into one graphics. This is again very slow due to Opacity and the sheer amount of line points

Show[
 Graphics[{Opacity[.01], Thickness[0.01], lines}],
 StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}],
 Frame -> True, FrameTicks -> True, PlotRange -> {{-2, 2}, {-2, 2}}, 
 PlotRangeClipping -> True
 ]

but the result looks somewhat nice and I believe it is what you are after

Please be aware that there is much room for improvement.

Line integral convolution

Another idea is to simply provide a test-function that tells you to which sink the starting point converges. This is close to the line function but doesn't create so many line-points.

test[{x0_?NumericQ, y0_?NumericQ}] :=
 Module[{x, y, t, res},
  res = NDSolveValue[{x'[t] == y[t], 
     y'[t] == -.15 y[t] + x[t] - x[t]^3, x[0] == x0, y[0] == y0}, {x, 
     y}, {t, 0, 100}];
  Norm[{-1, 0} - Through[res[100]]] < .01
  ]

Now, you could create a colored image that shows you the regions of convergence in different colors. We will split a color scheme into two halves and introduce some noise which will be an excellent starting point for a LIC

mat = ParallelTable[
   List @@ 
    ColorData["TemperatureMap", 
     Boole[test[{x, y}]]/2 + RandomReal[1/2]], {y, -2, 
    2, .01}, {x, -2, 2, .01}];

bg = LineIntegralConvolutionPlot[{{y, -.15 y + x - x^3}, 
   Image[Reverse@mat]}, {x, -2, 2}, {y, -2, 2}, RasterSize -> 512]

and then

Show[
 bg,
 StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}, 
  StreamStyle -> Gray]
 ]

and finally, note that you can use this test-function with RegionPlot or you can employ it to color your StreamPlot accordingly:

stP = StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}, 
  StreamColorFunction -> 
   Function[{x, y}, 
    If[test[{x, y}], ColorData["TemperatureMap", .9], 
     ColorData["TemperatureMap", .1]]], 
  StreamColorFunctionScaling -> False]

Show[bg, stP]

Speed

As mentioned in the comments, the real system has no analytical solution. The bottlenecks in the last LIC example are (1) the creation of the domain image and (2) the calculation of the LIC itself. As for (2) there is nothing much we can do. However, one can speed up (1) by implementing a simple numerical scheme that can be compiled. Here is a hack of a simple Runge-Kutta solver that can be called in parallel. If you don't have a C compiler then you need to remove the compile to C option:

f[{x_, y_}] := {y, -.15 y + x - x^3}
step[yn_, h_] := Module[{k1, k2, k3},
  k1 = f[yn];
  k2 = f[yn + h/2.0*k1];
  k3 = f[yn - h*k1 + 2.0 h*k2];
  yn + h (1./6 k1 + 4./6 k2 + 1./6 k3)
  ]

fc = Compile[{{x0, _Real, 0}, {y0, _Real, 0}},
    Module[{eps = 0.001, xn = x0, yn = y0, h = .1, maxSteps = 5000, 
      sink = -1, xtmp = 0., ytmp = 0.},
     Do[
      {xtmp, ytmp} = #;
      If[
       Norm[{xtmp, ytmp} - {xn, yn}] < eps,
       Break[]
       ];
      {xn, yn} = {xtmp, ytmp},
      maxSteps
      ];
     If[Norm[{-1, 0} - {xn, yn}] < .1, 1, 0]
     ], Parallelization -> True, RuntimeAttributes -> {Listable}, 
    CompilationTarget -> "C"
    ] &[FullSimplify[step[{xn, yn}, h]]]

For a region sampling of 512x512 points, the runtime of the call to NDSolve within the test function is about 100s, the same compution with fc only needs 4s.

mat = fc @@ Transpose[
     Table[{x, y}, {y, -2, 2, 4/511}, {x, -2, 2, 4/511}], {2, 3, 
      1}]; // AbsoluteTiming
Image[Reverse[mat]]

I promise it will be the last graphics, but I quite like the style. It's in remembrance of a good friend of mine, who had a fast line integral convolution for Mathematica long before Wolfram could even spell the word. He always liked flow fishes in vector fields (which are called "Drop" in Mathematica) and often combined them:

lic = LineIntegralConvolutionPlot[{{y, -.15 y + x - x^3}, 
   Colorize[ImageAdjust@Image[
      GaussianFilter[Reverse[mat + RandomReal[{-1, 1}, {512, 512}]], 
       4]]]},  {x, -2, 2}, {y, -2, 2}, RasterSize -> 512];
Show[
 lic,
 StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
  StreamStyle -> White, StreamMarkers -> "Drop"]
 ]

In Jumps in NDSolve results, I found the separatrices and used them to construct polygons to color the basins of attraction, but since it had little to do with the question asked, I left the code for it out. Below is a similar approach. The plot is just two polygons, clipped to the plot range.

eq = {y, -15/100 y + x - x^3} /. {v : x | y :> v[t]}  (* ode system *)
(*  {y[t], x[t] - x[t]^3 - (3 y[t])/20}  *)

equil = Solve[eq == 0]   (* equilibria *)
(*  {{x[t] -> -1, y[t] -> 0}, {x[t] -> 0, y[t] -> 0}, {x[t] -> 1, y[t] -> 0}}  *)

(* Picking the separatrices that are the boundaries of the basins takes some
 * inspection of the system: Reverse the velocity field (eq -> -eq), and the basins
 * should be the same; however, the eigenvalues below would change sign.
 *)
sepeq = Pick[N@equil, Det@D[eq, {{x[t], y[t]}}] /. equil // Sign, -1] (* unstable ones *)
sep0 = Flatten[  (* initial velocity vectors of separatrices *)
    {#, -#} &@ # & /@
     MapThread[
      Pick[#2, Sign[#1], -1] &,  (* picks eigenvectors heading toward equilibrium *)
      Eigensystem[#]
      ],
    1] & /@ (D[eq, {{x[t], y[t]}}] /. sepeq)
(*
  {{x[t] -> 0., y[t] -> 0.}}
  {{{-0.680151, 0.680151}, {0.733072, -0.733072}}}  
*)

ics = MapThread[ (* IC = initial velocity times a small number *)
   {x[t], y[t]} + 100 # 2 $MachineEpsilon & /@ #2 /. #1 &,
   {sepeq, sep0}
   ];

sepsols = Flatten[  (* separatrix integral curves *)
   Map[
    NDSolve[{{x'[t], y'[t]} == eq, {x[0], y[0]} == #,
       WhenEvent[Abs[x[t]] > 2.2 && Abs[y[t]] > 2.2, "StopIntegration"]},
      {x, y}, {t, 0, -100},   (* eigenvalue is neg. so integrate backwards *)
      "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}] &,
    ics,
    {2}
    ],
   2];

lines = Cases[  (* get points of separatrices *)
   Normal@ ParametricPlot[Evaluate[{x[t], y[t]} /. sepsols], {t, -100, 100}, 
     PlotPoints -> 100, PlotRange -> All],
   Line[p_] :> p, Infinity];

basins = Module[{i},  (* boundaries of basins of attraction *)
     i = First@Nearest[#[[2]] -> "Index", Last@#[[1]]];
     Join[#[[1]], #[[2, i ;; 1 ;; -1]]]
     ] & /@ {Reverse /@ lines, Reverse[Reverse /@ lines]};

Show[
   Graphics[{Opacity[0.5], Riffle[{Red, Green}, Polygon /@ basins]}],
   #,
   Options@#
   ] &@
     StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
      StreamPoints -> {{{{1, 1}, Darker@Red}, {{1, 1.6}, Darker@Green}, Automatic}},
      PlotRange -> 2.05]

---

Update:

We can use the other separatrix as one of the stream lines and add a few more colored appropriately.

sep2 = Flatten[{#, -#} &@# & /@
      MapThread[Pick[#2, Sign[#1], 1] &, Eigensystem[#]], 
     1] & /@
   (D[eq, {{x[t], y[t]}}] /. sepeq);
ics2 = MapThread[{x[t], y[t]} + 100 # 2 $MachineEpsilon & /@ #2 /. #1 &, {sepeq, sep2}];
sepsols2 = Flatten[
   Map[
    NDSolve[{{x'[t], y'[t]} == eq, {x[0], y[0]} == #, 
       WhenEvent[x[t]^2 + y[t]^2 > 0.25, "StopIntegration"]},
      {x, y}, {t, 0, 100}] &,
    ics2, {2}],
   2];

Show[Graphics[{Opacity[0.5], 
     Riffle[{Red, Green}, Polygon /@ basins]}], #, Options@#] &@
 StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
  StreamPoints -> {Join[
     Thread@{{x[t], y[t]} /.
        PadRight[sepsols2, {Automatic, 3}, 
         List /@ Thread[t -> (x["Domain"] /. sepsols2)[[All, 1, 2]]]],
       {Darker@Red, Darker@Green}},
     {{{0., -0.3}, Darker@Red}, {{0., -0.62}, Darker@Red},
      {{0., 0.3}, Darker@Green}, {{0., 0.62}, Darker@Green}, Automatic}
     ]},
  PlotRange -> 2.05]