Bifurcation diagrams for multiple equation systems

Question

I am interested in constructing a bifurcation diagram for some of my parameters (especially for β and γ) in the dynamical system given in the code below. I want to see how parameter changes affect the stability of the system, but I am having some trouble in estimating the critical points in Mathematica and then using them to get the critical point orbits. I found this post helpful,

Mathematica code for Bifurcation Diagram

but the post refers to a logistic-equation-like case. I do not how to apply it to my multiple equation case. Any help or pointers would be most appreciated

The following code incldes some wrong code intended to get the critical points, so my apologies in advance.

f = {L, \[Psi], d} /. NDSolve[{
  d'[t] == (1/
       z) ( \[Alpha] - \[Beta] (d[t]/Y[t]) - \[Gamma] L[
         t] - \[Delta] \[Psi][t] - \[Tau] Y[t]) + 
    0.75 (d[t]/Y[t]) d[t],
  L'[t] == (j (Y[t]/K) - e)/N,
  \[Psi]'[t] == p L[t] - 0.35,
  Y[t] == (I + 
      E + \[Alpha] - \[Beta] ( d[t]/Y[t]) - \[Gamma] L[
        t] - \[Delta] \[Psi][
        t])/(1 - (1 - \[Tau]) (\[Psi][
          t] + (1 - s) (1 - \[Psi][t]) - m z)),
  L[0] == 0.50, \[Psi][0] == 0.50, 
  d[0] == 0} /. { \[Alpha] -> 0.05, \[Gamma] -> 0.75, \[Delta] -> 
   0.75, X -> 1, \[Beta] -> 0.05, j -> 0.45, \[Tau] -> 0.35, 
  I -> 1, s -> 0.85, p -> 0.75, E -> 1, m -> 1.2, e -> 0.035, 
  z -> 5.5, K -> 3.5, N -> 4},
{L, \[Psi], d, Y}, {t, 0, 300}, MaxSteps -> 1000000 ][[1]]

l = {L[t], \[Psi][t], d[t]};

cps[L] = f /. Quiet[Solve[f[[L'[t]]] == 0], All]

Show[GraphicsArray[
Table[Plot[f[[i]][t], {t, 0, 300}, PlotRange -> All,
PlotStyle -> Blue, Filling -> 0,
AxesLabel -> TraditionalForm /@ {t, l[[i]]}, 
DisplayFunction -> Identity], {i, 3}]
]]
data2 = Table[Evaluate[f[[#]][t] & /@ {1, 2, 3}], {t, 0, 300}];
Export["data2.csv", data2]
ListLinePlot[Transpose@data2, Filling -> 0, 
AxesLabel -> {t, {L, \[Psi], d}}, PlotRange -> All]
Graphics3D[Point[data2], BoxRatios -> 1, 
AxesLabel -> TraditionalForm /@ l]
ParametricPlot3D[Evaluate[Append[#[t] & /@ f, Red]], {t, 0, 300}, 
BoxRatios -> {1, 1, 1}, PlotRange -> All, PlotPoints -> 1500, 
AxesLabel -> TraditionalForm /@ l]       

Answer 1

First, I have a few general comments. The simplest bifurcation diagrams for differential equations involve a single parameter in a single equation and there are several illustrations of these on the Wolfram Demonstrations site. More generally, of course, a bifurcation occurs when we see a qualitative change in the behavior of a system as some parameter changes. For a system of equations, we could use the eigensystem of the linearization of the system in the neighborhood of an equilibrium to identify bifuractions. You, however, appear to have four differential-algebraic equations with 16 parameters. Also, a number of these (K, I, E, N) are actually reserved symbols; I'd avoid that.

I don't think I'm going to wade into this system that I know nothing about but I can present some ideas to study this kind of thing in the context of an example that I understand, namely the Selkov model presented on the Demonstrations site: http://demonstrations.wolfram.com/HopfBifurcationInTheSelkovModel/

That demonstration shows you how to study the bifurcation using NDSolve. It might make sense to study it using the groovy new ParametricNDSolve. First, the system is the following:

$$x' = -x + a y + x^2 y$$ $$y' = b - a y - x^2 y$$

We can solve this as follows.

pf = ParametricNDSolveValue[{
  x'[t] == -x[t] + a*y[t] + x[t]^2 y[t],
  y'[t] == b - a*y[t] - x[t]^2*y[t],
  x[0] == 0, y[0] == 2},
 {x, y}, {t, 0, 100}, {a, b}];

We can investigate the behavior with respect to the parameters as follows. If you hold $a=0.1$ and let be range, a two Hopf bifurcations (changes from fixed point to cycle and back) should be evident.

Manipulate[
 Block[{$PerformanceGoal = "Quality"},
  Show[{
    ContourPlot[-x + a*y + x^2 y == 0, {x, 0, 3}, {y, 0, 3},
     ContourStyle -> Dashed],
    ContourPlot[b - a*y - x^2*y == 0, {x, 0, 3}, {y, 0, 3},
     ContourStyle -> Dashed],
    StreamPlot[{-x + a*y + x^2 y, b - a*y - x^2*y},
     {x, 0, 3}, {y, 0, 3},
     StreamStyle -> Directive[Opacity[0.5]]],
    ParametricPlot[
     Evaluate[Through[pf[a, b][t]]],
     {t, 0, 100}, PlotRange -> {{0, 3}, {0, 3}},
     PlotStyle -> Directive[Thickness[0.007], Black]]}]],
 {{a, 0.1}, 0, 1}, {b, 0, 1}]

Portrayed as an animation, this looks like so: