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, ψ, d} /.
NDSolve[{d'[
t] == (1/z) (α - β (d[t]/Y[t]) - γ L[t] - δ ψ[
t] - τ Y[t]) + 0.75 (d[t]/Y[t]) d[t],
L'[t] == (j (Y[t]/K) - e)/N, ψ'[t] == p L[t] - 0.35,
Y[t] == (I + E + α - β (d[t]/Y[t]) - γ L[t] - δ ψ[
t])/(1 - (1 - τ) (ψ[t] + (1 - s) (1 - ψ[t]) - m z)),
L[0] == 0.50, ψ[0] == 0.50,
d[0] == 0} /. {α -> 0.05, γ -> 0.75, δ -> 0.75,
X -> 1, β -> 0.05, j -> 0.45, τ -> 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, ψ, d, Y}, {t, 0, 300}, MaxSteps -> 1000000][[1]]
l = {L[t], ψ[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, ψ, 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]