Plotting phase plane in Filippov systems

Question

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

Answer 1

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}}]}]
]