# Bifurcation diagram of nonlinear system of ODE [duplicate]

Welcome to all, I try to apply most suggestion codes of a plot bifurcation diagram but I could not, I got an error. Is there anyone who can help me to understand how I can plot a bifurcation diagram for the following system:

r = 0.431201; β1 = 2.99*10^-6; β2 = 1;
α1 = 0.44257; α2 = 0.4 ; α3 = 2.99*10^-6;
k1 = 0.11; k2 = 0.99; c1 = 0.33; c2 = 0.66;
sys3 = {
N1'[t] == r N1[t] (1 - β1 N1[t]) - β2 N1[t] T[t] + c1 V[t] N1[t],
T'[t] == α1 T[t] (1 - α2 T[t]) + α3 N1[t] T[t] - c2 V[t] T[t],
V'[t] == k1 - k2 V[t]};
con3 = {N1 == 1, T == 1, V == 2};
s3 = NDSolve[{sys3, con3}, {N1, T, V}, {t, 0, 360}]


here, I want to plot the bifurction diagram of the system, also, the @zhk solve my problem for parametric solution code, but stil my qustion does not had a complete solution.

## marked as duplicate by zhk, MarcoB, bbgodfrey differential-equations StackExchange.ready(function() { if (StackExchange.options.isMobile) return; $('.dupe-hammer-message-hover:not(.hover-bound)').each(function() { var$hover = $(this).addClass('hover-bound'),$msg = $hover.siblings('.dupe-hammer-message');$hover.hover( function() { $hover.showInfoMessage('', { messageElement:$msg.clone().show(), transient: false, position: { my: 'bottom left', at: 'top center', offsetTop: -7 }, dismissable: false, relativeToBody: true }); }, function() { StackExchange.helpers.removeMessages(); } ); }); }); Mar 12 at 1:41

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

• What do you want on the x and y axes of this bifurcation diagram? Also, could you say what the model represents? – Chris K Mar 4 at 7:07
• I want to see the effect of equation 3 on the behaviour of N1and T this model describes the behaviour of normal and tumor cells. – sana alharbi Mar 4 at 9:41
• Thanks, but I'm not sure what you mean by "equation 3". Is there one particular parameter of interest? – Chris K Mar 4 at 9:45
• can study more than one parameter, actually, the behaviour of the model has affected by \beta_2, c_1, c_2 , \alhpha _3 and k2 but I do not if I can examine all these parameters where already I know that the values of them – sana alharbi Mar 4 at 9:51
• @sanaalharbi How about this? ParametricPlot3D[Evaluate@{N1[t], T[t], V[t]} /. s3, {t, 0, 360}, PlotRange -> {{0, 3}, {0, 3}, {0, 3}}] – zhk Mar 4 at 16:10

## 2 Answers

I think you are looking for something like this

r = 0.431201; \[Beta]1 = 2.99*10^-6;\[Alpha]1 = 0.44257; \[Alpha]2 = 0.4;

Manipulate[

eq1 = x'[t] == r x[t] (1 - \[Beta]1 x[t]) - \[Beta]2 x[t] y[t] + c1 z[t] x[t];
eq2 = y'[t] == \[Alpha]1 y[t] (1 - \[Alpha]2 y[t]) + \[Alpha]3 x[t] y[t] - c2 z[t] y[t];
eq3 = z'[t] == k1 - k2 z[t];

sol = NDSolveValue[{eq1, eq2, eq3, x == x0, y == y0, z == z0}, {x, y, z}, {t, 0, 360}];

ParametricPlot3D[{sol[][t], sol[][t], sol[][t]}, {t, 0, 360}, AxesLabel -> {x, y, z},
PlotRange -> All,  ImageSize -> {300, 300}],
{{c2, 0.4975, "c2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{c1, 0.2215, "c1"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{k1, 0.8677, "k1"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{k2, 0.9611, "k2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{\[Beta]2, 0.9817, "\[Beta]2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{\[Alpha]3 , 0.2291, "\[Alpha]3"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
Delimiter,
{{x0, 1,"x"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
{{y0, 1, "y"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
{{z0, 0, "z"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
ControlPlacement -> Left, ContinuousAction -> False] Note: Experiment with different values for the parameters and the initial conditions.

• THANKS. THAT IS GOOD FOR FIND A PARAMETRIC PLOT3D. WHAT ABOUT THE BIFRECTION AS SHOWN BY LINK@ZHK – sana alharbi Mar 7 at 4:58
• @sanaalharbi You already got an answer for that by ChrisK – zhk Mar 7 at 5:04

I don't have time for a full analysis of your problem, but here are some thoughts and a rough approach to a bifurcation diagram.

First, the $$V$$ equation is decoupled from the others, so you should be able to solve for its equilibrium value $$\hat V=k_1/k_2$$ and inject it into the other two equations. Then regrouping of terms would reduce the number of independent parameters and also allow phase-plane analysis.

Once you've done that, I think you'll find the normal ($$N_1$$) and tumor cells ($$T$$) follow a combination of the Lotka-Volterra predator-prey model with logistic growth of both prey and predators. I think you could show that if a positive equilibrium exists it is unique and asymptotically stable.

This justifies a rough approach to generating a bifurcation diagram, using FindRoot to find an equilibrium for a given parameter value, then looping over the parameter value using the previous answer as an initial guess.

N1res = Tres = Vres = {}; (* results to be accumulated here *)
{N1i, Ti, Vi} = {21000, 0.6, k1/k2}; (* first initial guess *)
Do[
eq = FindRoot[sys3 /. {_'[t] -> 0, var_[t] -> var},
(* remove [t] and set time derivatives equal to zero *)
{N1, N1i}, {T, Ti}, {V, Vi}];
(* store results *)
AppendTo[N1res, {k1, N1 /. eq}];
AppendTo[Tres, {k1, T /. eq}];
AppendTo[Vres, {k1, V /. eq}];
(* update initial guess *)
{N1i, Ti, Vi} = {N1, T, V} /. eq
, {k1, 0.6, 10.0, 0.1}]

GraphicsRow[{
ListPlot[N1res, AxesLabel -> {"k1", "N1"}],
ListPlot[Tres, AxesLabel -> {"k1", "T"}],
ListPlot[Vres, AxesLabel -> {"k1", "V"}]
}, ImageSize -> Large] This could be improved in many ways, but hopefully it gets you started.

• I analyzed the stability of the model it was stable when the cells succeed to live together. but there is no growth for the third equation.@chris k – sana alharbi Mar 4 at 11:26
• The equation for $V$ is independent of the other variables $N_1$ and $T$, so in all cases $V\to k_1/k_2$ as $t\to\infty$. – Chris K Mar 4 at 11:33
• yes, that means, V is bounded, I proved that anlytically – sana alharbi Mar 4 at 15:32