# 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] == 0.25,
x[0] == 0.35}, {x, y}, {t, 0, 0.1}] /. B01


How I can plot:

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

2. y[t] for y[0] 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[0] == y0,
x[0] == 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[97] /@ Range[11];

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[0], y[0]} == 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[[1]]@"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  }}],

Show[pl3d, PlotRange -> {0.02 {-1, 1}, All, All}]

• @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