I have four ODEs corresponding to four species:
dRdt = g*R (1 - R/k) - c1*M*R - c2*P*R;
dMdt = e1*c1*M*R - c3*Z*M*(1 - s) - M*d1;
dPdt = e2*c2*P*R - c5*Z*P*s - u1*P*(M + P) + u2*R*Z - P*d2;
dZdt = e3*c3*Z*M*(1 - s) + e5*c5*Z*P*s + u1*P*(M + P) - u2*R*Z - Z*d3;
I check that everything works by first simulating the dynamics over time:
sim = {R -> R[t], M -> M[t], P -> P[t], Z -> Z[t]};
dR = {R'[t] == First@{dRdt /. sim}, R[0] == 1};
dM = {M'[t] == First@{dMdt /. sim}, M[0] == 1};
dP = {P'[t] == First@{dPdt /. sim}, P[0] == 1};
dZ = {Z'[t] == First@{dZdt /. sim}, Z[0] == 1};
Assigning my parameters and solving with NDSolve:
pars = {g -> 5, k -> 25, c1 -> 1, c2 -> 0.5, c3 -> 0.5,
c5 -> 0.5, e1 -> 0.5, e2 -> 0.5, e3 -> 0.5, e5 -> 0.5, d1 -> 0.5,
d2 -> 0.5, d3 -> 0.5, u1 -> 0.5, u2 -> 0.5, s -> 0.7};
sol = NDSolve[{{dR, dM, dP, dZ} /. pars}, {R, M, P, Z}, {t, 0,
1000}];
And finally plotting the simulations across my two variables of interest:
Manipulate[
sol = NDSolve[{{dR, dM, dP, dZ} /. {g -> 5, k -> var1, c1 -> 0.5,
c2 -> 0.5, c3 -> 0.5, c5 -> 0.5, e1 -> 0.5, e2 -> 0.5,
e3 -> 0.5, e5 -> 0.5, d1 -> 0.5, d2 -> 0.5, d3 -> 0.5,
u1 -> 0.5, u2 -> 0.5, s -> var2}}, {R, M, P, Z}, {t, 0,
10000}];
Plot[Evaluate[{R[t], M[t], P[t], Z[t]} /.
First[sol]], {t, 10,
10000}, PlotStyle -> cols, PlotLegends -> Automatic,
PlotRange -> {{0, 1000}, {0, 10}}],
{{var1, 50}, 0, 100, 1},
{{var2, 0.7}, 0, 1, 0.01}
]
This is great because I can use the toggle to visualize when certain species go extinct vs. when they all coexist. Now, I want to show exactly that using RegionPlot. Where I'm stuck is that I'm not sure how to ask RegionPlot to evaluate the areas of all species having equilibrium greater than zero as being True. I've tried the following but I know this is wrong:
Table[
sol = NDSolve[{{dR, dM, dP, dZ} /. {g -> 5, k -> var1, c1 -> 0.5,
c2 -> 0.5, c3 -> 0.5, c5 -> 0.5, e1 -> 0.5, e2 -> 0.5,
e3 -> 0.5, e5 -> 0.5, d1 -> 0.5, d2 -> 0.5, d3 -> 0.5,
u1 -> 0.5, u2 -> 0.5, s -> var2}}, {R, M, P, Z}, {t, 0,
10000}];
RegionPlot[
Evaluate[{R[t], M[t], P[t], Z[t]} /.
sol] > 0, {0, 100, 1}, {0,
1, 0.01}],
{{var1, 50}, 0, 100, 1},
{{var2, 0.7}, 0, 1, 0.01}
]
Any help is greatly appreciated, thank you!