Mathematica code for bifurcation diagram in 3D

Question

Good day. I need help with the code in mathematica to plot the bifurcation diagram (e vs z*) or (e vs x*), for the system x'[t]=x(1-x)-k1xy,y'[t]=ey(1-y)+k2xy-k3yz,z'[t]=-k4z+ek3yz, where X*=(x*,y*,z*) is a fixed point of the system, with k4=0.1,k1=0.02,k2=0.1,k3=1.2, and ebelongs to [0,4].

Thank you

Answer 1

Looks like a little three-species food web model -- a perfect excuse to use my new EcoEvo package.

First, install the package (one-time only):

PacletInstall["https://github.com/cklausme/EcoEvo/releases/download/v1.0.1/EcoEvo-1.0.1.paclet"]

Then load the package and set the model for analysis with SetModel:

<< EcoEvo`

SetModel[{
  Pop[x] -> {Equation :> x[t] (1 - x[t]) - k1 x[t] y[t]},
  Pop[y] -> {Equation :> e y[t] (1 - y[t]) + k2 x[t] y[t] - k3 y[t] z[t]},
  Pop[z] -> {Equation :> -k4 z[t] + e k3 y[t] z[t]}
  }]

Set the parameters, then solve for the equilibria:

k4 = 0.1; k1 = 0.02; k2 = 0.1; k3 = 1.2;
eq = SolveEcoEq[]

Now you can plot the parts of the final, coexistence equilibrium vs e:

Plot[x /. eq[[-1]], {e, 0, 4}, PlotRange -> {0, All}]
Plot[z /. eq[[-1]], {e, 0, 4}, PlotRange -> {0, All}]

Just for fun, you can simulate the model with EcoSim:

e = 1;
sol = EcoSim[{x -> 1, y -> 1, z -> 0.1}, 100];
PlotDynamics[sol]