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]