# Plot with NDSolve for a range of initial values

I have a differential equation system as follow:

B01 = {P -> 70, R -> 5, W -> 200, E1 -> 500, C1 -> 300, S -> 1, H -> 800};

G11 = (E1 - R*P - C1)*S /. B01
G12 = (E1 - H)*S /. B01
G13 = -C1*S /. B01
G14 = -H*S /. B01

C11 = (R*P - W)*S /. B01
C12 = -W*S /. B01
C13 = 0;
C14 = 0;

UG1[t_] := y[t] (G11) + (1 - y[t]) (G13)
UG2[t_] := y[t] (G12) + (1 - y[t]) (G14)

UC1[t_] := x[t] (C11) + (1 - x[t]) (C12)
UC2[t_] := x[t] (C13) + (1 - x[t]) (C14)

Solution1 = NDSolve[{x'[t] == x[t] (1 - x[t]) (UG1[t] - UG2[t]), y'[t] ==  y[t] (1 - y[t]) (UC1[t] - UC2[t]), y == 0.25,
x == 0.35}, {x, y}, {t, 0, 0.1}] /. B01


How I can plot:

1. x[t] for x in range of [0, 1] with interval 0.1 (for example: x = {0, 0.1, 0.2, 03, ...})

2. y[t] for y in range [0, 1] with interval 0.1

3. and at the end a parametric plot of $$x$$ vs $$y$$ similar to the picture below: Clear["Global*"]

B01 = {P -> 70, R -> 5, W -> 200, E1 -> 500, C1 -> 300, S -> 1,
H -> 800};
G11 = (E1 - R*P - C1)*S;
G12 = (E1 - H)*S;
G13 = -C1*S;
G14 = -H*S;
C11 = (R*P - W)*S;
C12 = -W*S;
C13 = 0;
C14 = 0;
UG1[t_] := y[t] (G11) + (1 - y[t]) (G13)
UG2[t_] := y[t] (G12) + (1 - y[t]) (G14)
UC1[t_] := x[t] (C11) + (1 - x[t]) (C12)
UC2[t_] := x[t] (C13) + (1 - x[t]) (C14)

eqns = {x'[t] == x[t] (1 - x[t]) (UG1[t] - UG2[t]),
y'[t] == y[t] (1 - y[t]) (UC1[t] - UC2[t]), y == y0,
x == x0} /. B01 // Simplify;

Solution1 = ParametricNDSolve[eqns, {x, y},
{t, 0, 1/10}, {x0, y0}];

Plot[Evaluate[Flatten[Table[x[x0, y0][t],
{x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
{t, 0, 1/10},
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {t, x}),
AspectRatio -> 1] Plot[Evaluate[Flatten[Table[y[x0, y0][t],
{x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
{t, 0, 1/10},
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {t, y}),
AspectRatio -> 1] ParametricPlot[Evaluate[Flatten[Table[{x[x0, y0][t], y[x0, y0][t]},
{x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
{t, 0, 1/10},
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {x, y}),
AspectRatio -> 1,
PlotPoints -> 75,
MaxRecursion -> 5] EDIT: Coordinating the PlotStyles

styles = ColorData /@ Range;

Plot[
Evaluate[
Flatten[
Table[x[x0, y0][t],
{y0, 0, 1, 1/10}, {x0, 0, 1, 1/10}],
1] /. Solution1],
{t, 0, 1/10},
PlotStyle -> styles,
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {t, x}),
AspectRatio -> 1] Plot[
Evaluate[
Flatten[
Table[y[x0, y0][t],
{x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}],
1] /. Solution1],
{t, 0, 1/10},
PlotStyle -> styles,
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {t, y}),
AspectRatio -> 1] ParametricPlot[
Evaluate[
Flatten[
Table[{x[x0, y0][t], y[x0, y0][t]},
{x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}],
1] /. Solution1],
{t, 0, 1/10},
PlotStyle -> styles,
PlotRange -> All,
Frame -> True,
FrameLabel -> (Style[#, 14] & /@ {x, y}),
AspectRatio -> 1,
PlotPoints -> 75,
MaxRecursion -> 5] • Thanks Bob, This is what I need. Jul 7 at 5:29

A quick way to get #3 is using StreamPlot:

StreamPlot[{x[t] (1 - x[t]) (UG1[t] - UG2[t]),
y[t] (1 - y[t]) (UC1[t] - UC2[t])}, {x[t], 0, 1}, {y[t], 0, 1}] • It is amazing, thanks Chris. Jul 7 at 5:30

When I saw this question, I first thought StreamPlot, but it doesn't do the other graphs. Then I though wouldn't it be cool to convert an NDSolve solution to BezierCurve. I thought it would be cool because I had never done it. But after doing that, which might become a separate Q&A, I thought that the following would be cooler, even though it's not as robust ParametricNDSolveValue (see @BobHanlon's answer).

OP's setup in lowercase:

b01 = {p -> 70, r -> 5, w -> 200, e1 -> 500, c1 -> 300, s -> 1,
h -> 800};

g11 = (e1 - r*p - c1)*s /. b01;
g12 = (e1 - h)*s /. b01;
g13 = -c1*s /. b01;
g14 = -h*s /. b01;

c11 = (r*p - w)*s /. b01;
c12 = -w*s /. b01;
c13 = 0;
c14 = 0;

ug1[t_] := y[t] (g11) + (1 - y[t]) (g13);
ug2[t_] := y[t] (g12) + (1 - y[t]) (g14);

uc1[t_] := x[t] (c11) + (1 - x[t]) (c12);
uc2[t_] := x[t] (c13) + (1 - x[t]) (c14);

ics =(* we'll use the mesh coordinates as initial conditions *)
DiscretizeRegion[Disk[{1/2, 1/2}, 1/2],
MaxCellMeasure -> 0.01]; solution2 = NDSolveValue[
{x'[t] == x[t] (1 - x[t]) (ug1[t] - ug2[t]),
y'[t] == y[t] (1 - y[t]) (uc1[t] - uc2[t]),
(* we'll use the mesh coordinates all at once!: *)
{x, y} == Transpose@MeshCoordinates@ics},
{x, y}, {t, -0.07, 0.07}];
(* dimensions of Through[solution2@"ValuesOnGrid"]
* are {coord, step, ic} = 2 x 250 x 88 *)
systraj = Transpose[ (* transpose dims to {ic, step, coord} *)
Through[solution2@"ValuesOnGrid"], {3, 2, 1}] //
{Automatic, Automatic, 3}, (* {ic, step, new coords} *)
solution2[]@"Grid" (* use t from x sol, = same t as y *)
] &;
Dimensions@systraj
(*  {88, 250, 3}  *)


Graphics:

We can project the 3D trajectories onto 2D planes:

Table[
Graphics[{
Riffle[
colors,
Line /@ systraj[[All, All, parts]]]
}, Options@Plot],
{parts, {{1, 2}, {1, 3}, {2, 3}}}] Or we can project the trajectories onto planes in 3D:

proj // ClearAll;
proj[k_, offsets_ : {-0.1, -0.3, -0.3}] :=(*
project onto coordinate plane *)
TranslationTransform[
ReplacePart[{0., 0., 0.}, k -> offsets[[k]]]] .
ScalingTransform[ReplacePart[{1., 1., 1.}, k -> 0.]];

colors = Hue[ (* VertexColors for the ICs *)
Rescale[ArcTan @@ #, {-Pi, Pi}, {0, 1}],
2 Norm@#,
1
] &@(# - {1/2, 1/2}) & /@ MeshCoordinates@ics;

plall = Graphics3D[{
Opacity[0.5], Thickness@0.003,
Riffle[
colors,
Line /@ systraj],
Opacity[0.2],
Table[Riffle[
colors,
Line /@ proj[k]@systraj], {k, 3}],
Append[ (* Initial conditions Disk[] *)
ReplacePart[
#,
1 -> PadLeft[First@#, {Automatic, 3}, 0.]
] &@First@Show@ics /. {_Directive :> EdgeForm[Gray]},
VertexColors -> colors]
},
BoxRatios -> {1, 1, 1},
Axes -> True,
AxesLabel -> {t, x, y},
ViewPoint -> {2, 2.3, 1.8}, ViewVertical -> {0, 0, 1}] Like the plots in the OP:

pl3d = Graphics3D[{
Opacity[0.5], Thickness@0.003,
Riffle[
colors,
Line /@ systraj],
Opacity[0.2],
Append[(* Initial conditions Disk[] *)
ReplacePart[
#,
1 -> PadLeft[First@#, {Automatic, 3}, 0.]
] &@First@Show@ics /. {_Directive :> EdgeForm[Gray]},
VertexColors -> colors]
},
BoxRatios -> {1, 1, 1},
Axes -> True,
AxesLabel -> {t, x, y},
ViewPoint -> {2, 2.3, 1.8}, ViewVertical -> {0, 0, 1}];

Grid[
{{Graphics[Inset@#2, ImageSize -> 170],
Graphics[Inset@#1, ImageSize -> 340],
SpanFromLeft}, {Graphics[Inset@#3, ImageSize -> 170],
SpanFromAbove, SpanFromBoth}},
Frame -> All, Alignment -> Top, Spacings -> {0, 0}
] & @@
Table[Show[pl3d, opts],
{opts, Transpose@{
{All,  {{0, 0.07}, All, All},  {{0, 0.07}, All, All}}],
{ None,   {1, -1}, {1, -1}},
{{1, -1},  None,   {-1, 1}},
{{-1, 1}, {-1, 1} , None  }}], Without the projections, the graphics are smaller and easier to manipulate. If one zooms in on the initial condition disk, one can see the dynamics of the system better:

Show[pl3d, PlotRange -> {0.02 {-1, 1}, All, All}] • Whoa! I especially like that last one. Jul 7 at 23:52
• oh my God, I am very impressed with your solutions. Your grapghs are very interesting and amazing. I appreciate and thank you very much for sharing these answers. Jul 8 at 15:00
• I thank you again for your solutions and plots Michael, Actually in my model time is time and it cant be negative, so we cant have negative parts of time in plots. Jul 8 at 15:17
• @Ahmad Thanks & you're welcome. As for time: (1) Since your model is autonomous, time is arbitrary. Solutions are invariant under translation $t\mapsto t+t_0$. So any solution may be translated to a time interval that is nonnegative. (2) You can always change the interval back to {t, 0, 0.1}. Note in the second-to-last graphic, the x and y vs. t plots are for nonnegative t. Using a symmetric interval was the easiest way to make the x vs. y plot to show Chris K's StreamPlot`; using negative t does not invalidate the plot because the model is autonomous. (3)... Jul 8 at 15:58
• (3) My principal goal was to show some techniques that you or another user might adapt. And to make some pretty pictures. :) Jul 8 at 15:58