differential equations - Cleaner Approach to Plotting Multiple Solution Curves on Phase Portrait

I am trying to plot multiple solution curves onto a phase portrait. I used the approach shown here, but is there a cleaner approach?

For example, can we just add a parameter $n$, that randomly selects $n$ random initial conditions in each quadrant and plot those solutions curves? Of course, those random ICs would be in the domain of the solution.

  graph1 = StreamPlot[{4/3 x + 2/3 y, 1/3 x + 5/3 y}, {x, -8, 8}, {y, -8, 8}];

  gensol[x0_, y0_] := 

        NDSolve[{x'[t] == 4/3 x[t] + 2/3 y[t], y'[t] == 1/3 x[t] + 5/3 y[t],
          x[0] == x0, y[0] == y0}, {x, y}, {t, -8, 8}];

   sol[1] = gensol[4, 0];
   sol[2] = gensol[-1, 0];
   sol[3] = gensol[-3, 3];
   sol[4] = gensol[3, 3];
   sol[5] = gensol[3, -3];

   graph2 = Table[

    ParametricPlot[Evaluate[{x[t] /. sol[i][[1]], y[t] /. sol[i][[1]]}], {t,-8, 8}, 
      PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> Hue[i/5]], {i, 
       5}] // Flatten;

  Show[Join[graph2, {graph1}], ImageSize -> 200]

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$

Answer

It's not quite that simple because random initial conditions may not generate useful trajectories, i.e. trajectories that show different parts of the phase plane and don't lie too close to each other. I wrote a package to find suitable initial conditions, and plot those, which you can find here: https://github.com/cekdahl/PhasePortrait

Here's an example of what it generates for your system:

PhasePortrait[{
  x'[t] == 4/3 x[t] + 2/3 y[t],
  y'[t] == 1/3 x[t] + 5/3 y[t]
  }, {x, y}, t, {{-8, -8}, {8, 8}},
 PortraitDensity -> 5
 ]

PortraitDensity controls the number of trajectories. The higher the density, the more trajectories you get.

You can read more about the package in the read me file on Github.