# Plotting eigenvalues

I am trying to reproduce figure 10 in the linked paper. The authors have used Matlab bvp4c function to find multiple solutions for the system of ODE and then carried out the stability.

Here is my try,

s = 1; beta = 0.1; M = 0.2; Pr = 1; max = 7;

pfun = ParametricNDSolveValue[{f'''[x] + f[x]*f''[x] +
beta*(1 - f'[x]^2) + M^2*(1 - f'[x]) == 0,
1/Pr*T''[x] + f[x]*T'[x] == 0, f[0] == s, f'[0] == L1, T[0] == 1,
f'[xmax] == 1, T[xmax] == 0}, {f, T}, {x, 0, xmax}, {L1}]; (*Eqs 8,9 and 10*)

part[l_List, n_Integer] := Part[l, n];

qfun = ParametricNDSolveValue[{F'''[x] + part[pfun[L1], 1][2]*F''[x] +
D[part[pfun[L1], 1][2], x, x]*
F[x] - (2*beta*D[part[pfun[L1], 1][2], x] + M^2 - g1)*F'[x] ==
0, 1/Pr*G''[x] + part[pfun[L1], 1][2] + G'[x] +
D[part[pfun[L1], 2][2], x]*F[x] + g1*G[x] == 0, F[0] == 0,
F'[0] == 0, G[0] == 0, F'[xmax] == 0, G[xmax] == 0}, {F, G}, {x,
0, xmax}, {L1, g1}]; (*Eqs 22, 23 and 34 using solution from 8,9 and 10*)


Here the part function is calling the solutions from the original system of ODEs. The NDSolve is even unable to find the one solution, so I cannot proceed further. Any suggestions?

If you write

s = 1; beta = 0.1; M = 0.2; Pr = 1; xmax = 7; (* note xmax, not max *)


then

pfun =
ParametricNDSolveValue[
{f'''[x] + f[x]*f''[x] + beta*(1 - f'[x]^2) + M^2*(1 - f'[x]) == 0,
1/Pr*T''[x] + f[x]*T'[x] == 0,
f[0] == s, f'[0] == L1, T[0] == 1, f'[xmax] == 1, T[xmax] == 0},
{f, T}, {x, 0, xmax}, {L1}];


will return something that will plot:

With[{L1 = 10.},
Plot[Evaluate[Through[pfun[L1][x]]], {x, 0, xmax},
PlotRange -> {-1, All},
AxesOrigin -> {0, -1}]]


However, in the case of next system of ODE's, I find I need rewrite your code to

With[{L1 = 10.},
pf1 = pfun[L1][[1]]; pf2 = pfun[L1][[2]];
qfun =
ParametricNDSolveValue[{F'''[x] + pf1[2]*F''[x] +
D[pf1[x], x, x][[0]][2]*F[x] - (2*beta*D[pf1[x], x][[0]][2] +
M^2 - g1)*F'[x] == 0,
1/Pr*G''[x] + pf1[2] + G'[x] + D[pf2[x], x][[0]][2]*F[x] + g1*G[x] == 0,
F[0] == 0, F'[0] == 0, G[0] == 0, F'[xmax] == 0, G[xmax] == 0},
{F, G}, {x, 0, xmax}, {g1}]];


and then qfun can be plotted for L1 = 10. and g1 = 2 with this code:

With[{g1 = 2},
Plot[Evaluate[Through[qfun[g1][x]]], {x, 0, xmax},
PlotRange -> All,
AxesOrigin -> {0, -120}]]


### Update

Looking the conditions placed on F by OP equations and the resulting plot, I conclude that F is identically zero for all x in the domain {0, 7}. This allows a much simpler solution for G to be developed, from which is it rather easy to write a function that will plot g over the given domain for and given pair of the parameters L1 and g1. I present that function in this update.

GfunPlot[l1_, g1_] :=
Block[{pf1, Gfun},
pf1 = pfun[l1][[1]];
Gfun =
ParametricNDSolveValue[
{1/Pr*G''[x] + pf1[2] + G'[x] + g*G[x] == 0,
G[0] == 0, G[xmax] == 0},
G, {x, 0, xmax}, g];
Plot[Gfun[g1][x], {x, 0, xmax},
PlotRange -> All,
PlotLabel ->
Style[Row[{G[x], ":  L1 = ", l1, ", g1 = ", g1}], Black, 14]]]


And here is a sample set of plots.

GraphicsGrid[Table[GfunPlot[l1, g1], {l1, {5, 10}}, {g1, {2, 4}}],
ImageSize -> Large]


This is really all the information that can be extracted from the OP's formulation of the problem. The conditions imposed on F constricts the system's behavior severely.

• First of Thx. I am looking for the plot g1 vs L1, as given in the paper (see Fig 10 on page 7)
– zhk
Dec 8, 2019 at 8:17
• @zhk that looks to be essentially the first plot, at least from rough visual inspection. Dec 8, 2019 at 8:42
• @CATrevillian The results plotted here are for f, F, T and G vs x.
– zhk
Dec 8, 2019 at 11:06
• @zhk. But you formulated your question with g1 and L1 as parameters of systems set up to find f, T, F and G as functions of x. I can only work with what you post on this site. I am not going read your reference material and work out the math for you. I very likely wouldn't even understand the reference. I'm here to help people with Mathematica, no solve their math or physics problems. There are other SE sites for that. Dec 8, 2019 at 17:34