# Mathematica code for bifurcation diagram in 3D

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

• or the bifurcation diagram (e vs z), where z is a solution Sep 23, 2019 at 13:40
• Welcome to Mathematica.SE, Patrick! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. Sep 23, 2019 at 14:18

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] • There's a bug in your add-on: the documentation pages are asking for /Users/klaus/Library/Mathematica/SystemFiles/FrontEnd/StyleSheets/MyStyleSheet.nb` ! Sep 23, 2019 at 15:44
• @murray Ugh! I'll try to fix it tonight. Thanks for the note! Sep 23, 2019 at 16:43
• @ChrisK: I just copied what you posted and was able to duplicate the results. I did get one error message: "Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result." Thanks.
– Moo
Sep 23, 2019 at 21:34
• @moo That's more of a warning than an error. Sep 24, 2019 at 6:34
• @Moo I put a few examples in the built-in docs - check out EcoEvo/tutorial/ExampleModels in the help browser. Sep 24, 2019 at 12:13