# Bifurcation diagrams

I'm trying to get the Bifurcation diagrams for u, v, w, s, f, g and h as function of r

I got some kind of code to do this, but I do n't have any idea of what is going on with code and i'm getting bunch of errors

the problem is like:

Parameters

rf = 20000;r0 = 10.^-7;m = 2;k = 2;rE = 1000;a = 0.001;
a1 = (m (2 + m Sqrt[2 - n] - n +
E^((2 m)/Sqrt[2 - n]) (-2 + m Sqrt[2 - n] + n)) )/((1 + E^((2 m)/
Sqrt[2 - n])) (2 - n)^(3/2));
\[Tau] = 1/(\[Pi]^2 + k^2);
r = Ta;


The system of the algebraic equations

F1 = k^2 \[Tau] Ta v[r] - u[r]/\[Tau] + (n u[r] )/\[Tau] +
16/(3 \[Pi]) Ta k^2 \[Tau] v[r] s[r];
F2 = -(\[Pi]^2/n)  s[r] ;
F3 = -(\[Pi]^2/2) u[r] v[r] - (4 \[Pi]^2)/n w[r];
F4 = \[Pi] u[r] w[r] + u[r] - v[r]/(n \[Tau]) ;
F5 = rE/2 g[r] +
16/(6 \[Pi]) rE g[r] s[r] - (4 n)/a f[r] + (4 n^2)/a f[r];
F6 = - (rE/2) f[r] +
16/(6 \[Pi]) rE f[r] s[r] - (4 n)/a g[r] + (4 n^2)/a g[r];
F7 = - ((4 n)/a) h[r] + (4 n^2)/a h[r] - a1 rE u[r];


getting the maximum of each row

G1 = Table[{D[F1, r]}, {n, 0.2, 0.6, 0.1}]
Dimensions[G1]
G2 = Table[{D[F2, r]}, {n, 0.2, 0.6, 0.1}];
G3 = Table[{D[F3, r]}, {n, 0.2, 0.6, 0.1}];
G4 = Table[{D[F4, r]}, {n, 0.2, 0.6, 0.1}];
G5 = Table[{D[F5, r]}, {n, 0.2, 0.6, 0.1}];
G6 = Table[{D[F6, r]}, {n, 0.2, 0.6, 0.1}];
G7 = Table[{D[F7, r]}, {n, 0.2, 0.6, 0.1}];


Flatten

Difeq = Flatten[{Table[{Simplify[G1[[i, 1]]] == 0,
Simplify[G2[[i, 1]]] == 0, Simplify[G3[[i, 1]]] == 0,
Simplify[G4[[i, 1]]] == 0, Simplify[G5[[i, 1]]] == 0,
Simplify[G6[[i, 1]]] == 0, Simplify[G7[[i, 1]]] == 0}, {i, 1, 5,
1}]}, 1];
Dimensions[Difeq]

soln = Flatten[{Table[{Solve[{Difeq[[i, 1]], Difeq[[i, 2]],
Difeq[[i, 3]], Difeq[[i, 4]], Difeq[[i, 5]], Difeq[[i, 6]],
Difeq[[i, 7]]}, {Derivative[1][u][r], Derivative[1][v][r],
Derivative[1][w][r], Derivative[1][s][r], Derivative[1][f][r],
Derivative[1][g][r], Derivative[1][h][r]}]}, {i, 1, 5, 1}]},
3];
Dimensions[soln]
equation =
Flatten[{Table[{eq[i, 1] =
Simplify[Derivative[1][u][r] /. soln[[i, 1]]],
eq[i, 2] = Simplify[Derivative[1][v][r] /. soln[[i, 2]]],
eq[i, 3] = Simplify[Derivative[1][w][r] /. soln[[i, 3]]],
eq[i, 4] = Simplify[Derivative[1][s][r] /. soln[[i, 4]]],
eq[i, 5] = Simplify[Derivative[1][f][r] /. soln[[i, 5]]],
eq[i, 6] = Simplify[Derivative[1][g][r] /. soln[[i, 6]]],
eq[i, 7] = Simplify[Derivative[1][h][r] /. soln[[i, 7]]]}, {i,
1, 5, 1}]}, 1];
Dimensions[equation]
sol = Flatten[{Table[{NDSolve[{Derivative[1][u][r] ==
Simplify[Derivative[1][u][r] /. soln[[i, 1]]],
Derivative[1][v][r] ==
Simplify[Derivative[1][v][r] /. soln[[i, 2]]],
Derivative[1][w][r] ==
Simplify[Derivative[1][w][r] /. soln[[i, 3]]],
Derivative[1][s][r] ==
Simplify[Derivative[1][s][r] /. soln[[i, 4]]],
Derivative[1][f][r] ==
Simplify[Derivative[1][f][r] /. soln[[i, 5]]],
Derivative[1][g][r] ==
Simplify[Derivative[1][g][r] /. soln[[i, 6]]],
Derivative[1][h][r] ==
Simplify[Derivative[1][h][r] /. soln[[i, 7]]],
u[r0] == 0.0001, v[r0] == 0.0001, w[r0] == 0.0001,
s[r0] == 0.0001, f[r0] == 0.0001, g[r0] == 0.0001,
h[r0] == 0.0001}, {u[r], v[r], w[r], s[r], f[r], g[r],
h[r]}, {r, r0, rf},
Method -> {"Extrapolation",
Method -> "LinearlyImplicitEuler"}]}, {i, 1, 5, 1}]}, 3];
Dimensions[sol]


Plotting

q1 = Plot[Evaluate[ u[Ta] /. sol], {Ta, r0, rf},
AxesLabel -> {r,  Subscript[U, s][ r ]}, PlotLabel -> "(A)" ,
Frame -> True,
FrameLabel -> {Style["r", Bold, Blue],
Style["\!$$\*SubscriptBox[\(U$$, $$s$$]\)[ r ]", Bold, Blue]},
PlotStyle -> {Black, Blue, Red, Thick, Dashed}]

• There are a few minor issues to fix first. It's a good idea not to use capitals for your function and variable names. Your def of Difeq runs a table up to i=7, but your Gi only have 5 elements. Your def of soln, equation and sol have Tables indexed by {i, 1, 1, 1}, which is "i runs from 1 to 1 in steps of 1". Is this a placeholder for something else, or can you remove the Table entirely? The curly brackets outside your {Table[...]}s are superfluous. You don't appear to need the levelspec in Flatten, eg: in your def of Difeq, Flatten[stuff, 7] === Flatten[stuff]. Commented Aug 8, 2017 at 0:19
• @aardvark2012 thank u for ur comment. I have edited the code. fortunately, I don't get any errors now, but the plot that I receive is blank Commented Aug 8, 2017 at 10:57

A couple of changes can get the plot working. First, some expressions are undefined at r = 0 (try equation /. Ta -> 0). But that can easily be fixed by setting r0 = 10.^-7. Second, since you set r = Ta you need to replace r with Ta in your plot function:

q1 = Plot[Evaluate[u[Ta] /. sol], {Ta, r0, rf},
AxesLabel -> {r, Subscript[U, s][r]}, PlotLabel -> "(A)",
Frame -> True,
FrameLabel -> {Style["r", Bold, Blue],
Style["\!$$\*SubscriptBox[\(U$$, $$s$$]\)[ r ]", Bold, Blue]},
PlotStyle -> {Black, Blue, Red, Thick, Dashed}]


which gives

Hope that helps.