# Plotting phase plane in Filippov systems

I wanted to plot vector fields and nullclines for Modified Rosenzweig-MacArthur Model.
and I have a problem plotting p1 and p2 in the same graph on top of each other(I want the vector fields arrows in the range {x2,2,5,4}).
I also have another problem with plotting points of intersections between red and green curves.

a = 0.3556;
b = 0.33;
d = 0.0444;
e = 0.2067;

(*plotting vector field in S1 *)
p1 = StreamPlot[{x1*(1 - x1) - (a*x1*x2)/(b + x1), (a*x1*x2)/(
b + x1) - d*x2}, {x1, 0, 1}, {x2, 0, 2.5}, Axes -> True,
StreamStyle -> Gray, StreamColorFunction -> None];
(*plotting vector field in S2 *)
p2 = StreamPlot[{x1*(1 - x1) - (a*x1*x2)/(b + x1), (a*x1*x2)/(
b + x1) - (d + e)*x2}, {x1, 0, 1}, {x2, 2.5, 4}, Axes -> True,
StreamStyle -> Gray, StreamColorFunction -> None];

(*Solving equations and then plotting the solution trajectory *)
deq1 = x1'[t] == x1[t]*(1 - x1[t]) - (a*x1[t]*x2[t])/(b + x1[t]);
deq2 = x2'[t] == (a*x1[t]*x2[t])/(b + x1[t]) - d*x2[t];
solution1 =
NDSolve[{deq1, deq2, x1[0] == .0001, x2[0] == .0001}, {x1[t],
x2[t]}, {t, 0, 200}];
p3 = ParametricPlot[
Evaluate[{x1[t], x2[t]} /. solution1], {t, 0, 200},
PlotStyle -> Black, PlotStyle -> Thickness[0.005]];

(*plotting nullclines *)
f[x1_, x2_] = x1*(1 - x1) - (a*x1*x2)/(b + x1);
g1[x1_, x2_] = (a*x1*x2)/(b + x1) - d*x2;
g2[x1_, x2_] = (a*x1*x2)/(b + x1) - (d + e)*x2;

(*plotting nullclines of f1 *)
cp1 = ContourPlot[{f[x1, x2] == 0, g1[x1, x2] == 0}, {x1, 0, 1}, {x2,
0, 3}, ContourStyle -> {Green, Red}];
l1 = Graphics[{Thick, Green, Line[{{0, 0}, {0, 3}}]}];
(*plotting nullclines of f2 *)
cp2 = ContourPlot[{f[x1, x2] == 0, g2[x1, x2] == 0}, {x1, 0, 1}, {x2,
0, 3}, ContourStyle -> {Green, Red}];
l2 = Graphics[{Thick, Red, Line[{{0, 0}, {1, 0}}]}];

(*points of intersections *)
ptRules1 = NSolve[{f[x1, x2] == 0, g1[x1, x2] == 0}, {x1, x2}];
ptRules2 = NSolve[{f[x1, x2] == 0, g2[x1, x2] == 0}, {x1, x2}];

(*plotting sliding segment *)
l3 = Graphics[{Thick, Black, Line[{{0, 2.5}, {1, 2.5}}]}];

(*Showing the whole figure *)
Show[p1, p2, p3, l1, l2, l3, cp1, cp2, FrameLabel -> {"x1", "x2"},
Frame -> True]



I will appreciate any help
warm regards

• Did you forget the definition of a? It's undefined above, and that causes lots of errors. Oct 3, 2022 at 11:32
• Can you post the produced plot as well. What is the problem exactly? Oct 3, 2022 at 11:32
• @MichaelE2 am sorry I forgot to put it in the code here,I edited the code Oct 3, 2022 at 11:34
• @DrMrstheMonarch I will try to put an image for the plot Oct 3, 2022 at 11:35
• Does adding the PlotRange option produce what you want?: Show[p1, p2, p3, l1, l2, cp1, cp2, FrameLabel -> {"x1", "x2"}, Frame -> True, PlotRange -> All] Oct 3, 2022 at 11:38

This should be more straightforward with my EcoEvo package.

First, install package (once ever):

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]


Then, load package, set model and parameters:

<< EcoEvo

SetModel[{
Pop[x1] -> {Equation :> x1*(1 - x1) - (a*x1*x2)/(b + x1), Color -> Green},
Pop[x2] -> {Equation :> (a*x1*x2)/(b + x1) - d*x2, Color -> Red},
Parameters :> {a >= 0, b >= 0, d >= 0}
}]

a = 0.3556; b = 0.33; d = 0.0444;


Simulate and plot phase plane (S1):

sol1 = EcoSim[{x1 -> 0.0001, x2 -> 0.0001}, 200];

pp1 = Show[
PlotEcoPhasePlane[{x1, 0, 1}, {x2, 0, 2.5}],
RuleListPlot[FinalSlice[sol1, 100], PlotStyle -> {Black, Thickness[0.005]}]
]


Equilibria:

eq1 = FindEcoEq[{x1, 0, 1}, {x2, 0, 4}]
(* {{x1 -> 0, x2 -> 0}, {x1 -> 1., x2 -> 0}, {x1 -> 0.0470823, x2 -> 1.01048}} *)


Change d to d+e and do it again (S2):

d = 0.0444 + 0.2067;
pp2 = PlotEcoPhasePlane[{x1, 0, 1}, {x2, 2.5, 4}];
eq2 = FindEcoEq[{x1, 0, 1}, {x2, 0, 4}]
(* {{x1 -> 0, x2 -> 0}, {x1 -> 1., x2 -> 0}, {x1 -> 0.792947, x2 -> 0.65385}} *)


Finally, put them together:

l3 = Graphics[{Thick, Black, Line[{{0, 2.5}, {1, 2.5}}]}];
Show[pp1, pp2, l3, PlotRange -> {{0, 1}, {0, 4}}]


More elegant, just put the two systems together with If in the x2 equation.

SetModel[{
Pop[x1] -> {Equation :> x1*(1 - x1) - (a*x1*x2)/(b + x1), Color -> Green},
Pop[x2] -> {Equation :> (a*x1*x2)/(b + x1) - If[x2 < 2.5, d*x2, (d + e)*x2],
Color -> Red},
Parameters :> {a >= 0, b >= 0, d >= 0, e >= 0}
}]

a = 0.3556; b = 0.33; d = 0.0444; e = 0.2067;

sol = EcoSim[{x1 -> 0.0001, x2 -> 0.0001}, 200];

Show[
PlotEcoPhasePlane[{x1, 0, 1}, {x2, 0, 4}],
RuleListPlot[FinalSlice[sol, 100], PlotStyle -> {Black, Thickness[0.005]}],
Graphics[{Thick, Black, Line[{{0, 2.5}, {1, 2.5}}]}]
]
`

• Thanks a lot for your help Nov 9, 2022 at 14:09