# Shooting technique using NDSolve

Here I am attaching my code below. It took all the variable's values, but the graph window was empty. I am not able to recognize, where is the error exactly. Differential equations were also, checked so many times. If any further information is required, please feel free to ask. Can anyone please help me out?

Constant values

\[Epsilon]1 = 0.4;    \[Epsilon]2 = 0.8;   \[Rho]w = 1;
\[Mu]w = 1;    CPw = 1;    \[Kappa]w = 1;    \[Sigma]w = 1;
btw = 1;    bcw = 1;    \[Rho]k = 1;    \[Mu]k = 1;
CPk = 1;    \[Kappa]k = 1;    \[Sigma]k = 1;    \[Beta]tk = 1;
\[Beta]ck = 1;    t = 1;    L = 1;    \[Sigma] = 0.001;
R = 2.5;    H = 1;    Da1 = 0.5;    Da2 = 0.5;    Da3 = 0.5;    \[Beta] = 1;
M = 1;    Gr = 1;    Gc = 1;    \[Theta] = 0;
\[Delta] = 1;    Pr = 1;    Ec = 1;
a1 = 1;    a2 = 1.5;    a3 = 2;    f = 1;
Nb = 1;    Nt = 1;    Le = 1;    \[CapitalGamma] = 1;(*CHEMICAL REACTION PARAMETER*)
K1 = 0.001;     Bi = 1;    n = 1;


Corresponding Variables

\[CurlyTheta] = (\[Mu]w*\[Rho]k)/(\[Mu]k*\[Rho]w);   \[Sigma] = \[Sigma]k/\[Sigma]w;
\[Mu] = \[Mu]k/\[Mu]w;    \[Rho] = \[Rho]w/\[Rho]k;
\[Beta]T = \[Beta]tk/btw;    \[Beta]C = \[Beta]ck/bcw;
Cp = CPw/CPk;    \[Rho]c = (\[Rho]w*CPw)/(\[Rho]k*CPk);
\[Kappa] = \[Kappa]w/\[Kappa]k;    \[Alpha] = (\[Rho]w*CPw*\[Kappa]k)/(\[Kappa]w*\[Rho]k*CPk);    Db = 1; Subscript[D, t] = 1; K = 1; Pr1 = 21;
H1 = Sqrt[\[CurlyTheta]]*H;
R1 = \[CurlyTheta]*R;
M1 = Sqrt[\[Sigma]/\[Mu]]*M;
Gr1 = \[Beta]T*\[CurlyTheta]*Gr;
Gc1 = \[Beta]C*\[CurlyTheta]*Gc;
Ec1 = Cp*Ec;
Nt1 = Nt/(Subscript[D, t]*\[Alpha]);
Nb1 = Nb/(Db*\[Alpha]);
Le1 = \[Alpha]*Le*Db;
\[Gamma]1 = (\[CurlyTheta]*\[CapitalGamma])/K;
K2 = (\[CurlyTheta]*K1)/K;
D1 = (-(1/Da1) - M^2);
D2 = Gr*Sin[\[Theta]];
D3 = Gc*Sin[\[Theta]];
D4 = \[Delta]*(R/H)^2;
D5 = 1/Pr;
D6 = (-a1 + f)/Pr;
D7 = (Ec/Da1 + (M^2*Ec));
D8 = Nb/Pr;
D9 = Nt/Pr;
D10 = Le*Pr*H^2;


(* Some more constants *) All the constants are getting values properly, all of them I checked.

D11 = Le*Pr*R;
D12 = Nt/Nb;
D13 = \[CapitalGamma]*Le*Pr;
D14 = K1*Le*Pr;
D15 = \[Rho]*H1^2;
D16 = (-(1/Da2) - M1^2);
D17 = Gr1*Sin[\[Theta]];
D18 = Gc1*Sin[\[Theta]];
D19 = 1/Pr1;
D20 = \[Delta]*(R1/H1)^2;
D21 = ((-a2 + f)*\[Kappa])/Pr;
D22 = (Ec1/Da2 + Ec1*M1^2);
D23 = (\[Rho]*Cp*Nb1)/Pr1;
D24 = (\[Rho]*Cp*Nt1)/Pr1;
D25 = Le1*Pr1*H1^2;
D26 = Le1*Pr1*R1;
D27 = Nt1/Nb1;
D28 = \[Gamma]1*Le1*Pr1;
D29 = K2*Le1*Pr1;
D30 = (-(1/Da3) - M^2);
D31 = (-a3 + f)/Pr;
D32 = (Ec/Da3 + (M^2*Ec));


(* Ordinary Differential equations and boundary conditions *)

n1 =
NDSolve[{u10''[y] - (R*u10'[y]) + D1*u10[y] + D2*T10[y] +
D3*\[Phi]10[y] + H^2 == 0,
u11''[y] - (R*u11'[y]) + ((D1 - I*H^2)*u11[y]) + D2*T11[y] +
D3*\[Phi]11[y] + H^2 == 0,
((D5 - D4)*T10''[y]) - (R*T10'[y]) - (D6*
T10[y]) + (Ec*(u10'[y])^2) + (D7*(u10[y])^2) + (D8*\[Phi]10'[
y]*T10'[y]) + (D9*T10'[y]^2) == 0,
((D5 - D4)*T11''[y]) - (R*T11'[y]) - ((D6 + I*H^2)*T11[y]) + (2*
Ec*u11'[y]*u10'[y]) + (2*D7*u11[y]*
u10[y]) + (D8*(\[Phi]11'[y]*T10'[y] + \[Phi]10'[y]*
T11'[y])) + (2*D9*T11'[y]*T10'[y]) == 0,
\[Phi]10''[y] - (D11*\[Phi]10'[y]) - (D13*\[Phi]10[y]) +
D12*T10''[y] - (D14) == 0,
\[Phi]11''[
y] - (D11*\[Phi]11'[y]) - ((D13 + I*D10)*\[Phi]11[y]) +
D12*T11''[y] == 0,
((1 + 1/\[Beta])*u20''[y]) - (R1*u20'[y]) + (D16*u20[y]) + (D17*
T20[y]) + (D18*\[Phi]20[y]) + D15 == 0,
((1 + 1/\[Beta])*u21''[y]) - (R1*u21'[y]) + ((D16 - I*H1^2)*
u21[y]) + (D17*T21[y]) + (D18*\[Phi]21[y]) + D15 == 0,
((D19 - D20)*T20''[y]) - (R1*T20'[y]) - (D21*
T20[y]) + (Ec1*(1 + 1/\[Beta])*(u20'[y])^2) + (D22*(u20[
y])^2) + (D23*\[Phi]20'[y]*T20'[y]) + (D24*T20'[y]^2) ==
0,
((D19 - D20)*T21''[y]) - (R1*T21'[y]) - ((D21 + I*H1^2)*
T21[y]) + (2*Ec1*(1 + 1/\[Beta])*u20'[y]*u21'[y]) + (2*D22*
u20[y]*u21[
y]) + (D23*(\[Phi]20'[y]*T21'[y] + \[Phi]21'[y]*
T20'[y])) + (2*D24*T20'[y]*T21'[y]) == 0,
\[Phi]20''[y] - (D26*\[Phi]20'[y]) - (D28*\[Phi]20[y]) +
D27*T20''[y] - (D29) == 0,
\[Phi]21''[
y] - (D26*\[Phi]21'[y]) - ((D28 + I*D25)*\[Phi]21[y]) +
D27*T21''[y] == 0,
u30''[y] - (R*u30'[y]) + D30*u30[y] + D2*T30[y] +
D3*\[Phi]30[y] + H^2 == 0,
u31''[y] - (R*u31'[y]) + ((D30 - I*H^2)*u31[y]) + D2*T31[y] +
D3*\[Phi]31[y] + H^2 == 0,
((D5 - D4)*T30''[y]) - (R*T30'[y]) - (D31*
T30[y]) + (Ec*(u30'[y])^2) + (D32*(u30[
y])^2) + (D8*\[Phi]30'[y]*T30'[y]) + D9*T30'[y]^2 == 0,
((D5 - D4)*T31''[y]) - (R*T31'[y]) - ((D31 + I*H^2)*T31[y]) + (2*
Ec*u30'[y]*u31'[y]) + (2*D32*u30[y]*
u31[y]) + (D8*(\[Phi]30'[y]*T31'[y] + \[Phi]31'[y]*
T30'[y])) + (2*D9*T30'[y]*T31'[y]) == 0,
\[Phi]30''[y] - (D11*\[Phi]30'[y]) - (D13*\[Phi]30[y]) +
D12*T30''[y] - (D14) == 0,
\[Phi]31''[
y] - (D11*\[Phi]31'[y]) - ((D13 + I*D10)*\[Phi]31[y]) +
D12*T31''[y] == 0,
u10[0] == L*u10'[0], u11[0] == L*u11'[0], u30[1] == -L*u30'[1],
u31[1] == -L*u31'[1], u10[\[Epsilon]1] == u20[\[Epsilon]1],
u11[\[Epsilon]1] == u21[\[Epsilon]1],
u20[\[Epsilon]2] == u30[\[Epsilon]2],
u21[\[Epsilon]2] == u31[\[Epsilon]2],
u10'[\[Epsilon]1] == \[Mu]*u20'[\[Epsilon]1],
u11'[\[Epsilon]1] == \[Mu]*u21'[\[Epsilon]1],
u30'[\[Epsilon]2] == \[Mu]*u20'[\[Epsilon]2],
u31'[\[Epsilon]2] == \[Mu]*u21'[\[Epsilon]2],
T10'[0] == -Bi*(T10[0]),
T11'[0] == -Bi*(T11[0]), T30'[1] == -Bi*(T30[1] - 1),
T31'[1] == -Bi*(T31[1]), T10[\[Epsilon]1] == T20[\[Epsilon]1],
T11[\[Epsilon]1] == T21[\[Epsilon]1],
T20[\[Epsilon]2] == T30[\[Epsilon]2],
T21[\[Epsilon]2] ==
T31[\[Epsilon]2], \[Kappa]*T10'[\[Epsilon]1] ==
T20'[\[Epsilon]1], \[Kappa]*T11'[\[Epsilon]1] ==
T21'[\[Epsilon]1], \[Kappa]*T30'[\[Epsilon]2] ==
T20'[\[Epsilon]2], \[Kappa]*T31'[\[Epsilon]2] ==
T21'[\[Epsilon]2],
\[Phi]10'[0] == -n*(\[Phi]10[0]),
\[Phi]11'[0] == -n*(\[Phi]11[0]), \[Phi]30'[
1] == -n*(\[Phi]30[1] - 1), \[Phi]31'[
1] == -n*(\[Phi]31[
1]), \[Phi]10[\[Epsilon]1] == \[Phi]20[\[Epsilon]1], \
\[Phi]11[\[Epsilon]1] == \[Phi]21[\[Epsilon]1], \[Phi]20[\[Epsilon]2] \
== \[Phi]30[\[Epsilon]2], \[Phi]21[\[Epsilon]2] == \
\[Phi]31[\[Epsilon]2],
Db*\[Phi]10'[\[Epsilon]1] == \[Phi]20'[\[Epsilon]1],
Db*\[Phi]11'[\[Epsilon]1] == \[Phi]21'[\[Epsilon]1], \[Phi]20'[\
\[Epsilon]2] == Db*\[Phi]30'[\[Epsilon]2], \[Phi]21'[\[Epsilon]2] ==
Db*\[Phi]31'[\[Epsilon]2]}, {u10, u20, u30, u11, u21, u31, T10,
T20, T30, T11, T21,
T31, \[Phi]10, \[Phi]20, \[Phi]30, \[Phi]11, \[Phi]21, \
\[Phi]31}, {y, 0, 1},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton",
"StepControl" -> "LineSearch"}}];


(* Plot command *)

p1 = Plot[
Piecewise[{{Re[(u30[y] + (\[Sigma]*(E^(I*t))*u31[y])) /.
n1], \[Epsilon]2 <= y <=
1}, {Re[(u20[y] + (\[Sigma]*(E^(I*t))*u21[y])) /.
n1], \[Epsilon]1 <=
y <= \[Epsilon]2}, {Re[(u10[
y] + (\[Sigma]*(E^(I*t))*u11[y])) /. n1],
0 <= y <= \[Epsilon]1}}], {y, 0, 1}, PlotRange -> All]

• Your data uploaded without semicolon. Please, put ; in the end in every position with definitions like D32 = (Ec/Da3 + (M^2*Ec)); Jan 31, 2023 at 10:05
• @AlexTrounev done sir Jan 31, 2023 at 10:16
• In few minutes we have a message NDSolve::nderr: Error test failure at y == 0.8557735490281189; unable to continue. With options PrecisionGoal -> 5, AccuracyGoal -> 5 code generates some solution, but also with errors. Jan 31, 2023 at 13:05
• Yes sir, is there any way sir we can rectify this problem @AlexTrounev Jan 31, 2023 at 15:00
• Several of the ODEs can be solved symbolically with DSolve. Do that first to reduce the dimensionality of the problem. Feb 1, 2023 at 1:06

Original solution

With a bit of assistance, NDSolve can solve these ODEs using the Shooting method. To begin, set the List of ODEs and boundary conditions in the question to eqs and then split it into three parts.

eq0 = Cases[eqs, z_ /; MemberQ[z, T10[_] | T20[_] | T30[_] | ϕ10[_] |
ϕ20[_] | ϕ30[_] | T10'[_] | T20'[_] | T30'[_] |
ϕ10'[_] | ϕ20'[_] | ϕ30'[_], Infinity] &&
! MemberQ[z, T11[_] | T21[_] | T31[_] | ϕ11[_] | ϕ21[_] |
ϕ31[_] | T11'[_] | T21'[_] | T31'[_] | ϕ11'[_] |
ϕ21'[_] | ϕ31'[_], Infinity]];

eq1 = Cases[eqs, z_ /; MemberQ[z, T11[_] | T21[_] | T31[_] | ϕ11[_] |
ϕ21[_] | ϕ31[_] | T11'[_] | T21'[_] | T31'[_] |
ϕ11'[_] | ϕ21'[_] | ϕ31'[_], Infinity]];

equ = Complement[eqs, eq0, eq1];


The point of doing so is that equ can be solved symbolically.

su = DSolve[equ, {u10, u20, u30, u11, u21, u31}, {y, 0, 1}] //
Simplify // Flatten;


Next, use su to eliminate {u10, u20, u30, u11, u21, u31} from eq0 and solve it numerically to obtain solutions for {T10, T20, T30, ϕ10, ϕ20, ϕ30}.

Chop@Simplify[eq0 /. su];
s0 = NDSolve[%, {T10, T20, T30, ϕ10, ϕ20, ϕ30}, {y, 0, 1},
Method -> {"Shooting", "StartingInitialConditions" ->
{T10[ϵ2] == 0, T20[ϵ2] == 0, ϕ10[ϵ2] == 0, ϕ20[ϵ2] == 0,
T10'[ϵ2] == 0, T20'[ϵ2] == 0, ϕ10'[ϵ2] == 0, ϕ20'[ϵ2] == 0}}] // Flatten;


Note that the Shooting starts from y = ϵ2 instead of the default y = 0. to avoid the error message given in the question. Finally, eliminate all the variables obtained so far to solve eq1 numerically.

Chop@Simplify[eq1 /. su /. s0];
s1 = NDSolve[%, {T11, T21, T31, ϕ11, ϕ21, ϕ31}, {y, 0, 1},
Method -> {"Shooting", "StartingInitialConditions" ->
{T11[ϵ2] == 0, T21[ϵ2] == 0, ϕ11[ϵ2] == 0, ϕ21[ϵ2] == 0,
T11'[ϵ2] == 0, T21'[ϵ2] == 0, ϕ11'[ϵ2] == 0, ϕ21'[ϵ2] == 0}}] // Flatten;


Plots of the solutions are given by

ReImPlot[Evaluate@Through[su[[All, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


Plot[Evaluate@Through[s0[[All, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


ReImPlot[Evaluate@Through[s1[[All, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


That the solutions satisfy the boundary conditions can be verified, if desired, by back-substitution. For instance,

(Subtract @@@ eq0[[7 ;;]]) /. s0
(* {1.54431*10^-9, -2.29796*10^-11, 5.62184*10^-10, 0.,
1.27355*10^-8, 0., 2.1361*10^-8, -2.97741*10^-10, 2.46916*10^-9, 0.,
7.96958*10^-8, 0.} *)


Second, simpler solution

sall = NDSolve[eqs, {u10, u20, u30, u11, u21, u31, T10, T20,
T30, ϕ10, ϕ20, ϕ30, T11, T21, T31, ϕ11, ϕ21, ϕ}, {y, 0, 1},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" ->
{T10[ϵ2] == 0, T20[ϵ2] == 0, ϕ10[ϵ2] == 0, ϕ20[ϵ2] == 0, T10'[ϵ2] == 0,
T20'[ϵ2] == 0, ϕ10'[ϵ2] == 0, ϕ20'[ϵ2] == 0, T11[ϵ2] == 0, T21[ϵ2] == 0,
ϕ11[ϵ2] == 0, ϕ21[ϵ2] == 0, T11'[ϵ2] == 0, T21'[ϵ2] == 0, ϕ11'[ϵ2] == 0,
ϕ21'[ϵ2] == 0, u10[ϵ2] == 0, u20[ϵ2] == 0, u10'[ϵ2] == 0, u20'[ϵ2] == 0,
u11[ϵ2] == 0, u21[ϵ2] == 0, u11'[ϵ2] == 0, u21'[ϵ2] == 0}}] // Flatten


Then,

ReImPlot[Evaluate@Through[sall[[;; 6, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
ReImPlot[Evaluate@Through[sall[[7 ;; 12, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
ReImPlot[Evaluate@Through[sall[[13 ;;, 2]][y]], {y, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


reproduces the plots above. So, why does this work while the similar approach in the question does not? A possibility is that shooting from y = ϵ2 instead of the default y = 0 requires fewer elements (24 instead of 30) for "StartingInitialConditions". Also, y = ϵ2 is closer to where a few of the solutions become large (of order 350). A better response, however, might be, "Who knows?" I have found the error, "nderr: Error test failure" to occur somewhat randomly in NDSolve calculations. The only suggestion in the documentation, to try Method -> "StiffnessSwitching" does not help here.

• Thank you so much, @bbgodfre sir, for the help, but this is very difficult for me to understand. Can you please explain it in easy language? Feb 1, 2023 at 9:49
• @KomalGoyal The second solution, which I added a few hours ago, probably is easier to follow Feb 2, 2023 at 0:00
• @bbgodfrey Second solution is very good (+1). Feb 7, 2023 at 14:38
• @bbgodfrey, Please can you intervene in the post bellow for the numerical resolution of a difficult algebraic equation, mathematica.stackexchange.com/q/279693/58233 Feb 11, 2023 at 0:03