# Solution of non linear coupled equations by using shooting, block and StartingInitialConditions

My code is not giving out put. Moreover i do not know how to copy paste this code from my note book. is there any easiest way to paste here. I m plotting graph for the different values of R. I have also to plot for the differnt values of epsilon......... can we predict values given in block. ............................................................................................................................................................................................................................................................................................................................................................................................................................................................ kindly guide

eqn1=(1+1/\[Beta])((1+\[Epsilon]-\[Epsilon] \[Theta][y]) f''''[y]-\[Epsilon] f''[y] \[Theta]''[y]-2 \[Epsilon] f'''[y] \[Theta]'[y])+\[Alpha] (y f'''[y]+3 f''[y])+R (f[y] f'''[y]-f'[y] f''[y])- M^2 f''[y]==0;
eqn2=(1+NN)\[Theta]''[y]+Pr  (\[Alpha](m \[Theta][y]+y \[Theta]'[y])+ R(f[y] \[Theta]'[y]-m f'[y] \[Theta][y]))==0;
bcs1={f[-1]==A,f'[-1]==0, f[1]==1, f'[1]==0 };
bcs2={\[Theta][-1]==1, \[Theta][1]==0};
f1=Block [{\[Beta]=0.5, \[Epsilon]=1.5, Pr=1, M=1,NN=2,R=-5,m=1, \[Alpha]=0.5,A=-0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2},
{f[\[Eta]],f'[\[Eta]], f''[\[Eta]],f'''[\[Eta]], \[Theta][\[Eta]],\
[Theta]'[\[Eta]]}, {\[Eta],-1,1},
Method->{"Shooting",  "StartingInitialConditions" ->
{{f[-1]==A,f[-1]==0,f'[-1]==0, f''[-1]==0,f'''[-1]==0,\[Theta][-1]==1,\
[Theta][-1]==0,\[Theta]'[-1]==0}}}
]]
];
TableForm[
Table[{\[Eta],{f[\[Eta]]/.f1}, {f'[\[Eta]]/.f1}, {f''[\[Eta]]/.f1},{\[Theta]
[\[Eta]]/.f1},{\[Theta]'[\[Eta]]/.f1}},{\[Eta],-1,1,0.2}], TableAlignments->
{Center,Center},
[Theta][\[Eta]]","\[Theta]'[\[Eta]]"}},TableSpacing->{1,5}]
f2=Block [{\[Beta]=0.5, \[Epsilon]=1.5, Pr=1, M=1,NN=2,R=-2.5,m=1, \
[Alpha]=0.5,A=-0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2},
{f[\[Eta]],f'[\[Eta]], f''[\[Eta]],f'''[\[Eta]], \[Theta][\[Eta]],\
[Theta]'[\[Eta]]}, {\[Eta],-1,1},
Method->{"Shooting",  "StartingInitialConditions" ->
{{f[-1]==A,f[-1]==0,f'[-1]==0, f''[-1]==0,f'''[-1]==0,\[Theta][-1]==1,\
[Theta][-1]==0,\[Theta]'[-1]==0}}}
]]
];
TableForm[
Table[{\[Eta],{f[\[Eta]]/.f2}, {f'[\[Eta]]/.f2}, {f''[\[Eta]]/.f2},{\
[Theta][\[Eta]]/.f2},{\[Theta]'[\[Eta]]/.f2}},{\[Eta],-1,1,0.2}],
TableAlignments->{Center,Center},
[Theta][\[Eta]]","\[Theta]'[\[Eta]]"}},TableSpacing->{1,5}]
f3=Block [{\[Beta]=0.5, \[Epsilon]=1.5, Pr=1, M=1,NN=2,R=0,m=1, \
[Alpha]=0.5,A=-0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2},
{f[\[Eta]],f'[\[Eta]], f''[\[Eta]],f'''[\[Eta]], \[Theta][\[Eta]],\
[Theta]'[\[Eta]]}, {\[Eta],-1,1},
Method->{"Shooting",  "StartingInitialConditions" ->
{{f[-1]==A,f[-1]==0,f'[-1]==0, f''[-1]==0,f'''[-1]==0,\[Theta][-1]==1,\
[Theta][-1]==0,\[Theta]'[-1]==0}}}
]]
];
TableForm[
Table[{\[Eta],{f[\[Eta]]/.f3}, {f'[\[Eta]]/.f3}, {f''[\[Eta]]/.f3},{\
[Theta]
[\[Eta]]/.f3},{\[Theta]'[\[Eta]]/.f3}},{\[Eta],-1,1,0.2}], TableAlignments->
{Center,Center},
[Theta][\[Eta]]","\[Theta]'[\[Eta]]"}},TableSpacing->{1,5}]
f4=Block [{\[Beta]=0.5, \[Epsilon]=1.5, Pr=1, M=1,NN=2,R=2.5,m=1, \
[Alpha]=0.5,A=-0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2},
{f[\[Eta]],f'[\[Eta]], f''[\[Eta]],f'''[\[Eta]], \[Theta][\[Eta]],\
[Theta]'[\[Eta]]}, {\[Eta],-1,1},
Method->{"Shooting",  "StartingInitialConditions" ->
{{f[-1]==A,f[-1]==0,f'[-1]==0, f''[-1]==0,f'''[-1]==0,\[Theta][-1]==1,\
[Theta][-1]==0,\[Theta]'[-1]==0}}}
]]
];
TableForm[
Table[{\[Eta],{f[\[Eta]]/.f4}, {f'[\[Eta]]/.f4}, {f''[\[Eta]]/.f4},{\[Theta]
[\[Eta]]/.f4},{\[Theta]'[\[Eta]]/.f4}},{\[Eta],-1,1,0.2}], TableAlignments->
{Center,Center},
[Theta][\[Eta]]","\[Theta]'[\[Eta]]"}},TableSpacing->{1,5}]
f5=Block [{\[Beta]=0.5, \[Epsilon]=1.5, Pr=1, M=1,NN=2,R=5,m=1, \
[Alpha]=0.5,A=-0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2},
{f[\[Eta]],f'[\[Eta]], f''[\[Eta]],f'''[\[Eta]], \[Theta][\[Eta]],\
[Theta]'[\[Eta]]}, {\[Eta],-1,1},
Method->{"Shooting",  "StartingInitialConditions" ->
{{f[-1]==A,f[-1]==0,f'[-1]==0, f''[-1]==0,f'''[-1]==0,\[Theta][-1]==1,\
[Theta][-1]==0,\[Theta]'[-1]==0}}}
]]
];
TableForm[
Table[{\[Eta],{f[\[Eta]]/.f5}, {f'[\[Eta]]/.f5}, {f''[\[Eta]]/.f5},{\[Theta]
[\[Eta]]/.f5},{\[Theta]'[\[Eta]]/.f5}},{\[Eta],-1,1,0.2}], TableAlignments->
{Center,Center},
[Theta][\[Eta]]","\[Theta]'[\[Eta]]"}},TableSpacing->{1,5}]
Needs["PlotLegends"]
Plot[Evaluate[{{f[\[Eta]]}/.f1,f[\[Eta]]/.f2,f[\[Eta]]/.f3,f[\[Eta]]/.f4,f[\
[Eta]]/.f5}],
{\[Eta],-1,1},ImageSize->500,PlotRange->All,
PlotStyle->{{Dotted,Black},{DotDashed,Red},{Dashed,Green},{Blue}},Axes->
{False,False},
Frame->True,FrameLabel->{Style["\[Eta]",Italic,Black],Style["f(\
[Eta])",Italic,Black]},
PlotLegend->{"R=-5","R=-2.5", "R=0","R=2.5","R=5"},LegendPosition->
{-0.1,-0.0},
Plot[Evaluate[{{f'[\[Eta]]}/.f1,f'[\[Eta]]/.f2,f'[\[Eta]]/.f3,f'[\
[Eta]]/.f4,f'[\[Eta]]/.f5}],
{\[Eta],-1,1},ImageSize->500,PlotRange->All,
PlotStyle->{{Dotted,Black},{DotDashed,Black},{Dashed,Brown},{Gray}},Axes->
{False,False},
Frame->True,FrameLabel->{Style["\[Eta]",Italic,Black],Style["f'(\
[Eta])",Italic,Black]},
PlotLegend->{"R=-5","R=-2.5", "R=0","R=2.5","R=5"},LegendPosition->
{-0.1,-0.0},
Plot[Evaluate[{{\[Theta][\[Eta]]}/.f1,\[Theta][\[Eta]]/.f2,\[Theta][\
[Eta]]/.f3,\[Theta][\[Eta]]/.f4,\[Theta][\[Eta]]/.f5}],
{\[Eta],-1,1},ImageSize->500,PlotRange->All,
PlotStyle->{{Dotted,Black},{DotDashed,Black},{Dashed,Brown},{Gray}},Axes->
{False,False},
Frame->True,FrameLabel->{Style["\[Eta]",Italic,Black],Style["\[Theta](\
[Eta])",Italic,Black]},
PlotLegend->{"R=-5","R=-2.5", "R=0","R=2.5","R=5"},LegendPosition->
{-0.1,-0.0},


The above code is giving message that my dependent varibles are more than equations. here is the copied message NDSolve::underdet

• Equations eqn1 and eqn2 depend on y, and you call the solver with the $\eta$ variable. Fix the whole code. Replace $\eta$ in NDSolve[] with y. Sep 2, 2019 at 15:26

In version 12, a solution can be obtained without specifying a method. It is only necessary to replace $$\eta$$ with y in NDSolve.

eqn1 = (1 +
1/\[Beta]) ((1 + \[Epsilon] - \[Epsilon] \[Theta][y]) f''''[
y] - \[Epsilon] f''[y] \[Theta]''[y] -
2 \[Epsilon] f'''[y] \[Theta]'[y]) + \[Alpha] (y f'''[y] +
3 f''[y]) + R (f[y] f'''[y] - f'[y] f''[y]) - M^2 f''[y] == 0;
eqn2 = (1 + NN) \[Theta]''[y] +
Pr (\[Alpha] (m \[Theta][y] + y \[Theta]'[y]) +
R (f[y] \[Theta]'[y] - m f'[y] \[Theta][y])) == 0;
bcs1 = {f[-1] == A, f'[-1] == 0, f[1] == 1, f'[1] == 0};
bcs2 = {\[Theta][-1] == 1, \[Theta][1] == 0};
f1 = Block[{\[Beta] = 0.5, \[Epsilon] = 1.5, Pr = 1, M = 1, NN = 2,
R = -5, m = 1, \[Alpha] = 0.5, A = -0.2},
First[NDSolve[{eqn1, eqn2, bcs1, bcs2}, {f[y], f'[y], f''[y],
f'''[y], \[Theta][y], \[Theta]'[y]}, {y, -1, 1}]]]
Plot[Evaluate[{f[y] /. f1, f'[y] /. f1, \[Theta][y] /. f1}], {y, -1,
1}, ImageSize -> 500, PlotRange -> All,
PlotStyle -> {{Dotted, Black}, {DotDashed, Red}, {Dashed,
Green}, {Blue}}, Axes -> {False, False}, Frame -> True,
FrameLabel -> {Style["\[Eta]", Italic, Black],
Style["f(\[Eta]), f'(\[Eta]), \[Theta](\[Eta])", Italic, Black]},
PlotLegends -> Automatic]
`

• dear thanks allot. your given code is in mathematica 12 and its giving output at epsilon=1.5 but i am facing another problem while i have changed eta by y. code is not not running for the higher values of epsilon in mathematica 11. Kindly help in this version Sep 2, 2019 at 18:57
• Sorry, I only have versions 8, 9, 10 and 12. I'll see what can be done with growth epsilon. Is this a heat transfer problem? Sep 3, 2019 at 3:40
• yes its a heat transfer problem Sep 4, 2019 at 19:07