# How to Make a Bifurcation plot of a non linear system of ODEs

I am modelling bacterial growth. I have a quite difficult non linear systems of ODE's. I have 2 bacteria that i'm modelling with competition, as wel as transfer to a different place. so there are 2 ODEs for every bacteria because when it moves it will be a different one (P en S equation). Now i am trying to do a parameter analysis and i've been trying for days to make a bifurcation plot with on the x-axis one of the parameters (a0,b0,d1,k,awmin,tmin). But nothing i've been trying so far worked. Can you guys help me out a little and show how it works? here is the code of the bacteria, the bifurcation should take a spot around t=150, I've had a lot of trouble with this.

sp1 = D[G[t], t] ==
Growth*G[t]*((Sp - G[t] - a0*B[t] )/Sp) + P - S  ;
sp2 = D[M[t],
t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t] )/Ss) + S -
P   ;
P := M[t]* d1 ;
S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k     ;
sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t] )/Sp) + P2 - S2;
sp4 = D[A[t], t] ==
Growth2* A[t]*((Ss - A[t] - b0*M[t] )/Ss) + S2 - P2;
P2 := A[t]* d1;
S2 := B[t]*0.01*((B[t] + b0*G[t] )/Sp) + B[t]*k;
(*Parameters*)

Growth = 14.8*((T - tmin)*(Aw - Amin))^2  ;
Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2  ;
T := 4 + 18.09*Sin[0.01016*t + 0.3418];
Aw := 0.97 + 0.0001*t   ;
Kill = 1.5;
tmin = 12;
Amin = 0.95;
Sp = 1000;
Ss = 500;
a0 = 0.4;
b0 = 0.21;
d1 = 0.05;
k = 0.01;
s = 0.001;
all = { sp1, sp2, sp3, sp4};
init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30};
Solution =
NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}];
Plot[Solution, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}]
ParametricPlot[{Solution[[1]], Solution[[3]]}, {t, 0, 200},
PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red,
AspectRatio -> 1, PlotLabel -> "Phase plot", AxesLabel -> {"G", "B"}]


• This doesn't fix your problem, but you'd be better off if you rewrote your code in a more ordered fashion like this pastebin.com/LWVXcmWC where definitions come before usage and things are segregated into functions of t vs. parameters. Aug 21, 2020 at 17:42

We can use Module for parametric research as follows

f[a0p_, b0p_, d1p_, kp_, Ap_, tp_] :=
Module[{a0 = a0p, b0 = b0p, d1 = d1p, k = kp, Amin = Ap, tmin = tp},
Kill = 1.5;
(*tmin=12;
Amin=0.95;*)
Sp = 1000;
Ss = 500;
(*a0=0.4;
b0=0.21;
d1=0.05;
k=0.01;*)
s = 0.001;
sp1 = D[G[t], t] == Growth*G[t]*((Sp - G[t] - a0*B[t])/Sp) + P - S;
sp2 = D[M[t],
t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t])/Ss) + S - P;
P := M[t]*d1;
S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k;
sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t])/Sp) + P2 - S2;
sp4 = D[A[t], t] ==
Growth2*A[t]*((Ss - A[t] - b0*M[t])/Ss) + S2 - P2;
P2 := A[t]*d1;
S2 := B[t]*0.01*((B[t] + b0*G[t])/Sp) + B[t]*k;
Growth = 14.8*((T - tmin)*(Aw - Amin))^2;
Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2;
T := 4 + 18.09*Sin[0.01016*t + 0.3418];
Aw := 0.97 + 0.0001*t; all = {sp1, sp2, sp3, sp4};
init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30};
Solution =
NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}];
Solution]


To plot Solution we use f as

sol = f[.4, .21, .05, .01, .95, 12];
{Plot[sol, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}],
ParametricPlot[{sol[[1]], sol[[3]]}, {t, 0, 200},
PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red,
AspectRatio -> 1, PlotLabel -> "Phase plot",
AxesLabel -> {"G", "B"}]}


To plot Solution as a function of parameters for a given time t=150 we define new function, for example,

solp[ap_] := f[ap, .21, .05, .01, .95, 12] /. t -> 150;
solp1[bp_] := f[.4, bp, .05, .01, .95, 12] /. t -> 150;
var = {G, M, B, A}; Table[
Plot[solp[ap][[i]], {ap, .1, .5}, Frame -> True,
FrameLabel -> {"a0", var[[i]]}], {i, 4}]

Table[Plot[solp1[bp][[i]], {bp, .1, .5}, Frame -> True,
FrameLabel -> {"b0", var[[i]]}], {i, 4}]


• I think ParametricNDSolve(Value)[] might be more conducive to exploring parameter spaces. Aug 23, 2020 at 14:50
• @J.M. Yes, we can use ParametricNDSolve. It gives actually the same result. Aug 23, 2020 at 15:06
• Thank you so much for your response, it is very useful! You just made my day :D
– Zara
Aug 24, 2020 at 12:31
• @Zara You are welcome! Aug 24, 2020 at 14:09