# Solving a system of ODE using shooting technique with NDSolve command

In this code, I am trying to solve a system of first-order ODE with corresponding boundary conditions. But it is showing some error which I am not able to rectify. Can anyone help me in rectifying the errors of this code. I will be very thankful to you. If any doubt is there please feel free to ask.

Constants

$Beta] = 0.001; \[Epsilon] = 0.1; \[Alpha] = 1; \[Zeta] = 0; H = 1; M = 1; k = 1.5; z = 0.5*\[Pi]; Subscript[\[Phi], 1] = 0.01; Subscript[\[Phi], 2] = 0.01; \[Rho]w = 997.1; \[Mu]w = 0.001004; \[Sigma]w = 5.5*10^-6; \[Rho]k = 783; \[Mu]k = 0.00151; \[Sigma]k = 5*10^-11; \[Rho]fe = 5200; \[Sigma]fe = 25000; \[Rho]Mi = 5060; \[Sigma]Mi = 2.09*10^4; t = \[Pi]/3; Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^2.5; Subscript[A, 2] = (1 - Subscript[\[Phi], 1]) + (( Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]w); Subscript[A, 3] = ((\[Sigma]fe*(1 + (2*Subscript[\[Phi], 1]))) + (2*\[Sigma]fe*(1 - Subscript[\[Phi], 1])))/((\[Sigma]fe*(1 - Subscript[\[Phi], 1])) + (\[Sigma]w*(2 + Subscript[\[Phi], 1]))); Subscript[B, 1] = 1/((1 - Subscript[\[Phi], 2])^2.5*(1 - Subscript[\[Phi], 1])^2.5); Subscript[B, 2] = (((1 - Subscript[\[Phi], 1]) + (( Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 - Subscript[\[Phi], 2])) + ( Subscript[\[Phi], 2]*\[Rho]Mi)/\[Rho]k; Subscript[B, 3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w* Subscript[A, 3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 - Subscript[\[Phi], 2])) + (\[Sigma]w*Subscript[A, 3]*(2 + Subscript[\[Phi], 2]))); \[Sigma] = \[Sigma]w/\[Sigma]k; \[Mu] = \[Mu]w/\[Mu]k; \[Rho] = \[Rho]w/\[Rho]k; M1 = Sqrt[\[Mu]/\[Sigma]]*M; H1 = Sqrt[\[Mu]/\[Sigma]]*H;  Using NDSolve commands I am entering all the system of first order differential equations with corresponding boundary conditions n1 = NDSolve[{(Subscript[A, 1]/Subscript[A, 2]*(u100''[ y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2*u100[y]))) - (( Subscript[A, 3]*M^2)/Subscript[A, 2]*u100[y]) - (k/(y + k)* H^2/Subscript[A, 2]) == 0, (Subscript[A, 1]/Subscript[A, 2]*(u110''[ y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2*u110[y]))) - ((( Subscript[A, 3]*M^2)/Subscript[A, 2] + I*H^2)* u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 0, (Subscript[B, 1]/Subscript[B, 2]*(u200''[ y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2*u200[y]))) - (( Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]*u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 0, (Subscript[B, 1]/Subscript[B, 2]*(u210''[ y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2*u210[y]))) - ((( Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2] + I*H1^2)*u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[ B, 2]) == 0, (Subscript[A, 1]/Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[ y]*((Subscript[A, 1]/Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[ y]*((Subscript[A, 1]/Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[ y]*((Subscript[A, 1]/Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[ y]*((Subscript[A, 1]/Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[B, 1]/Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[ y]*((Subscript[B, 1]/Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + (( Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[ y]*((Subscript[B, 1]/Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + (( Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[ y]*((Subscript[B, 1]/Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + (( Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[ B, 2] + (I*H1^2)))) == 0, (Subscript[B, 1]/Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[ y]*((Subscript[B, 1]/Subscript[B, [![2$*(1/(y + k)^2 + $Alpha$^2)) + ((
Subscript$B, 3$*Subscript$A, 3$*$Sigma$*M1^2)/Subscript$B, 2$ + (I*H1^2)))) ==
0, (Subscript$A, 1$/Subscript$A, 2$*(g1''$y$ + (1/(y + k)*g1'$y$))) - (g1$y$*((Subscript$A, 1$/Subscript$A, 2$*1/(y + k)^2) + ((][1]][1]
Subscript[A, 3]*M^2)/Subscript[A, 2]))) ==
0, (Subscript[A, 1]/Subscript[A,
2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[
y]*((Subscript[A, 1]/Subscript[A,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((
Subscript[A, 3]*M^2)/Subscript[A, 2]))) ==
0, (Subscript[A, 1]/Subscript[A,
2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[
y]*((Subscript[A, 1]/Subscript[A,
2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2]))) ==
0, (Subscript[A, 1]/Subscript[A,
2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[
y]*((Subscript[A, 1]/Subscript[A, 2]*1/(y + k)^2) + ((
Subscript[A, 3]*M^2)/Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[A, 1]/Subscript[A,
2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[
y]*((Subscript[A, 1]/Subscript[A,
2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[A, 1]/Subscript[A,
2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[
y]*((Subscript[A, 1]/Subscript[A,
2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[
y]*((Subscript[B, 1]/Subscript[B, 2]*(1/(y + k)^2)) + ((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B,
2]))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[
y]*((Subscript[B, 1]/Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B,
2]))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[
y]*((Subscript[B, 1]/Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B,
2]))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[
y]*((Subscript[B, 1]/Subscript[B,

2]*1/(y + k)^2) + (((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[
B, 2]) + (I*H1^2)))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[
y]*((Subscript[B, 1]/Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[
B, 2]) + (I*H1^2)))) ==
0, (Subscript[B, 1]/Subscript[B,
2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[
y]*((Subscript[B, 1]/Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((
Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[
B, 2]) + (I*H1^2)))) == 0, u100[1] == 0, u110[1] == 0,
u200[-1] == 0, u210[-1] == 0, f2[1] == 0, f4[1] == 0,
f1[1] == -u100'[1], f3[1] == -u110'[1],
f5[-1] == -(Cos[\[Zeta]]*u200'[-1]),
f6[-1] == -(Sin[\[Zeta]]*u200'[-1]),
f7[-1] == -(Cos[\[Zeta]]*u210'[-1]),
f8[-1] == -(Sin[\[Zeta]]*u210'[-1]),
g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]),
g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]),
g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]),
g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]),
g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4*
u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]),
g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]),
g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4*
u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]),
g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), u100[0] == u200[0],
u110[0] ==
u210[0], ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(u100'[0] -u100[0]/k) == (u200'[0] -u200[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(u110'[0] - u110[0]/k) == (u210'[0] -u210[0]/k),
f1[0] == f5[0], f2[0] == f6[0], f3[0] == f7[0],
f4[0] ==
f8[0], ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(f1'[0] - f1[0]/k) == (f5'[0] - f5[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(f2'[0] - f2[0]/k) == (f6'[0] - f6[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(f3'[0] - f3[0]/k) == (f7'[0] - f7[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(f4'[0] - f4[0]/k) == (f8'[0] - f8[0]/k), g1[0] == g7[0],
g2[0] == g8[0], g3[0] == g9[0], g4[0] == g10[0], g5[0] == g11[0],
g6[0] ==
g12[0], ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g1'[0] - g1[0]/k) == (g7'[0] - g7[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g2'[0] - g2[0]/k) == (g8'[0] - g8[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g3'[0] - g3[0]/k) == (g9'[0] - g9[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g4'[0] - g4[0]/k) == (g10'[0] - g10[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g5'[0] - g5[0]/k) == (g11'[0] - g11[0]/
k), ((\[Mu]*Subscript[A, 1])/Subscript[B,
1])*(g6'[0] - g6[0]/k) == (g12'[0] - g12[0]/k)}, {u100, u110,
u200, u210, f1, f2, f3, f4, f5, f6, f7, f8, g1, g2, g3, g4, g5,
g6, g7, g8, g9, g10, g11, g12}, {y, -1, 1},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton",
"StepControl" -> "LineSearch"}} ];


Expressions

u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1;
u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1;
u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1;
u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1;
u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] +
g3[y]*Cos[2*\[Alpha]*z]) /. n1;
u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] +
g6[y]*Cos[2*\[Alpha]*z]) /. n1;
u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] +
g9[y]*Cos[2*\[Alpha]*z]) /. n1;
u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] +
g12[y]*Cos[2*\[Alpha]*z]) /. n1;
u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. n1;
u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. n1;
u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. n1;
u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. n1;


Plot Command for the expression

p1 = Plot[
Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <=
0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))],
0 <= y <= 1}}], {y, -1, 1}, PlotRange -> All, AspectRatio -> 1,
AxesLabel -> {y, U}, PlotStyle -> {Green, Thick}, ImageSize -> 260,
LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}]

• Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematica Meta, or in Mathematica Chat. Comments continuing discussion may be removed.
– Kuba
Jun 13, 2023 at 6:26

To solve this problem we need to express boundary conditions as follows

\[Beta] = 0.001;
\[Epsilon] = 0.1;
\[Alpha] = 1;
\[Zeta] = 0;
H = 1;
M = 1;
k = 1.5;
z = 0.5*\[Pi];
Subscript[\[Phi], 1] = 0.01;
Subscript[\[Phi], 2] = 0.01;
\[Rho]w = 997.1;
\[Mu]w = 0.001004;
\[Sigma]w = 5.5*10^-6;
\[Rho]k = 783;
\[Mu]k = 0.00151;
\[Sigma]k = 5*10^-11;
\[Rho]fe = 5200;
\[Sigma]fe = 25000;
\[Rho]Mi = 5060;
\[Sigma]Mi = 2.09*10^4;
t = \[Pi]/3;
Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^2.5;
Subscript[A,
2] = (1 -
Subscript[\[Phi], 1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]w);
Subscript[A,
3] = ((\[Sigma]fe*(1 + (2*
Subscript[\[Phi], 1]))) + (2*\[Sigma]fe*(1 -
Subscript[\[Phi], 1])))/((\[Sigma]fe*(1 -
Subscript[\[Phi], 1])) + (\[Sigma]w*(2 +
Subscript[\[Phi], 1])));
Subscript[B, 1] =
1/((1 - Subscript[\[Phi], 2])^2.5*(1 - Subscript[\[Phi], 1])^2.5);
Subscript[B,
2] = (((1 -
Subscript[\[Phi],
1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 -
Subscript[\[Phi], 2])) + (Subscript[\[Phi],
2]*\[Rho]Mi)/\[Rho]k;
Subscript[B,
3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w*
Subscript[A,
3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 -
Subscript[\[Phi], 2])) + (\[Sigma]w*
Subscript[A, 3]*(2 + Subscript[\[Phi], 2])));
\[Sigma] = \[Sigma]w/\[Sigma]k;
\[Mu] = \[Mu]w/\[Mu]k;
\[Rho] = \[Rho]w/\[Rho]k;
M1 = Sqrt[\[Mu]/\[Sigma]]*M;
H1 = Sqrt[\[Mu]/\[Sigma]]*H;

eq = {(Subscript[A, 1]/
Subscript[A,
2]*(u100''[
y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2*
u100[y]))) - ((Subscript[A, 3]*M^2)/Subscript[A, 2]*
u100[y]) - (k/(y + k)*H^2/Subscript[A, 2]) ==
0, (Subscript[A, 1]/
Subscript[A,
2]*(u110''[
y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2*
u110[y]))) - (((Subscript[A, 3]*M^2)/Subscript[A, 2] +
I*H^2)*u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) ==
0, (Subscript[B, 1]/
Subscript[B,
2]*(u200''[
y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2*
u200[y]))) - ((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
M1^2)/Subscript[B, 2]*
u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) ==
0, (Subscript[B, 1]/
Subscript[B,
2]*(u210''[
y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2*
u210[y]))) - (((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
M1^2)/Subscript[B, 2] + I*H1^2)*
u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[
y]*((Subscript[A, 1]/
Subscript[A,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2]))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[
y]*((Subscript[A, 1]/
Subscript[A,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2]))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[
y]*((Subscript[A, 1]/
Subscript[A,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[
y]*((Subscript[A, 1]/
Subscript[A,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/
Subscript[B, 2] + (I*H1^2)))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/
Subscript[B, 2] + (I*H1^2)))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g1''[y] + (1/(y + k)*g1'[y]))) - (g1[
y]*((Subscript[A, 1]/
Subscript[A, 2]*1/(y + k)^2) + (Subscript[A, 3]*M^2)/
Subscript[A, 2])) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[
y]*((Subscript[A, 1]/
Subscript[A,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[A, 3]*
M^2)/Subscript[A, 2]))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[
y]*((Subscript[A, 1]/
Subscript[A, 2]*(1/(y + k)^2 +
4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2]))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[
y]*((Subscript[A, 1]/
Subscript[A, 2]*1/(y + k)^2) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[
y]*((Subscript[A, 1]/
Subscript[A, 2]*(1/(y + k)^2 +
4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[A, 1]/
Subscript[A, 2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[
y]*((Subscript[A, 1]/
Subscript[A, 2]*(1/(y + k)^2 +
4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
Subscript[A, 2] + (I*H^2)))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[
y]*((Subscript[B, 1]/
Subscript[B, 2]*(1/(y + k)^2)) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[
y]*((Subscript[B, 1]/
Subscript[B,
2]*1/(y + k)^2) + (((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
H1^2)))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
H1^2)))) ==
0, (Subscript[B, 1]/
Subscript[B, 2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[
y]*((Subscript[B, 1]/
Subscript[B,
2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
H1^2)))) == 0};

bc = {u100[1] == 0, u110[1] == 0, u200[-1] == 0, u210[-1] == 0,
f2[1] == 0, f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1],
f5[-1] == -(Cos[\[Zeta]]*u200'[-1]),
f6[-1] == -(Sin[\[Zeta]]*u200'[-1]),
f7[-1] == -(Cos[\[Zeta]]*u210'[-1]),
f8[-1] == -(Sin[\[Zeta]]*u210'[-1]),
g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]),
g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]),
g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]),
g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]),
g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4*
u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]),
g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]),
g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4*
u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]),
g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), u100[0] == u200[0],
u110[0] ==
u210[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(u100'[0] -
u100[0]/k) == (u200'[0] -
u200[0]/k), ((\[Mu]*Subscript[A, 1])/
Subscript[B, 1])*(u110'[0] - u110[0]/k) == (u210'[0] -
u210[0]/k), f1[0] == f5[0], f2[0] == f6[0], f3[0] == f7[0],
f4[0] ==
f8[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f1'[0] -
f1[0]/k) == (f5'[0] -
f5[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f2'[0] -
f2[0]/k) == (f6'[0] -
f6[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f3'[0] -
f3[0]/k) == (f7'[0] -
f7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f4'[0] -
f4[0]/k) == (f8'[0] - f8[0]/k), g1[0] == g7[0], g2[0] == g8[0],
g3[0] == g9[0], g4[0] == g10[0], g5[0] == g11[0],
g6[0] ==
g12[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g1'[0] -
g1[0]/k) == (g7'[0] -
g7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g2'[0] -
g2[0]/k) == (g8'[0] -
g8[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g3'[0] -
g3[0]/k) == (g9'[0] -
g9[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g4'[0] -
g4[0]/k) == (g10'[0] -
g10[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g5'[0] -
g5[0]/k) == (g11'[0] -
g11[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g6'[0] -
g6[0]/k) == (g12'[0] - g12[0]/k)};

var = {u100, u110, u200, u210, f1, f2, f3, f4, f5, f6, f7, f8, g1, g2,
g3, g4, g5, g6, g7, g8, g9, g10, g11, g12}; coord = {y, -1,
1}; op =
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}};

var0 = Table[var[[i]]@y, {i, Length[var]}]; var1 = D[var0, y]; var2 =
D[var1, y];

sol0 = Solve[eq, var2];

sol01 = sol0 /. y -> 1; sol02 = sol0 /. y -> -1;

bc1 = bc /. sol01[[1]] /. sol02[[1]];


Solution

sol = NDSolve[{eq, bc1}, var, coord, op];


Visualization

Plot[Evaluate[Re[var0] /. sol[[1]]], coord, PlotRange -> All,
PlotLegends -> var, Frame -> True, FrameLabel -> Automatic]


Expressions definition

With[{n1 = sol[[1]]},
u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1;
u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1;
u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1;
u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1;
u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] +
g3[y]*Cos[2*\[Alpha]*z]) /. n1;
u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] +
g6[y]*Cos[2*\[Alpha]*z]) /. n1;
u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] +
g9[y]*Cos[2*\[Alpha]*z]) /. n1;
u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] +
g12[y]*Cos[2*\[Alpha]*z]) /. n1;
u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /.
n1;
u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /.
n1;
u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /.
n1;
u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /.
n1];


Plot U

p1 = Plot[
Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <=
0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))],
0 <= y <= 1}}], {y, -1, 1}, PlotRange -> All, AspectRatio -> 1,
AxesLabel -> {"y", "U"}, PlotStyle -> {Green, Thick},
ImageSize -> 260,
LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}]


Update 1. In a case of homogeneous boundary conditions for U we define first U[-1], U[1] as follows

U[-1] = ((u200[
y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] +
f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] +
g8[y]*Sin[2*\[Alpha]*z] +
g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] +
f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] +
g11[y]*Sin[2*\[Alpha]*z] +
g12[y]*Cos[2*\[Alpha]*z])))) /. y -> -1;

U[1] = ((u100[
y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] +
f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] +
g2[y]*Sin[2*\[Alpha]*z] +
g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] +
f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] +
g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /.
y -> 1;


If we suppose that U[-1]=U[1]=0 then we have new bc, for example

bc = {U[1] == 0, u110[1] == 0, U[-1] == 0, u210[-1] == 0, f2[1] == 0,
f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1],
f5[-1] == -(Cos[\[Zeta]]*u200'[-1]),
f6[-1] == -(Sin[\[Zeta]]*u200'[-1]),
f7[-1] == -(Cos[\[Zeta]]*u210'[-1]),
f8[-1] == -(Sin[\[Zeta]]*u210'[-1]),
g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]),
g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]),
g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]),
g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]),
g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4*
u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]),
g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4*
u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]),
g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]),
g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4*
u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]),
g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4*
u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), u100[0] == u200[0],
u110[0] ==
u210[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(u100'[0] -
u100[0]/k) == (u200'[0] -
u200[0]/k), ((\[Mu]*Subscript[A, 1])/
Subscript[B, 1])*(u110'[0] - u110[0]/k) == (u210'[0] -
u210[0]/k), f1[0] == f5[0], f2[0] == f6[0], f3[0] == f7[0],
f4[0] ==
f8[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f1'[0] -
f1[0]/k) == (f5'[0] -
f5[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f2'[0] -
f2[0]/k) == (f6'[0] -
f6[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f3'[0] -
f3[0]/k) == (f7'[0] -
f7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f4'[0] -
f4[0]/k) == (f8'[0] - f8[0]/k), g1[0] == g7[0], g2[0] == g8[0],
g3[0] == g9[0], g4[0] == g10[0], g5[0] == g11[0],
g6[0] ==
g12[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g1'[0] -
g1[0]/k) == (g7'[0] -
g7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g2'[0] -
g2[0]/k) == (g8'[0] -
g8[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g3'[0] -
g3[0]/k) == (g9'[0] -
g9[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g4'[0] -
g4[0]/k) == (g10'[0] -
g10[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g5'[0] -
g5[0]/k) == (g11'[0] -
g11[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g6'[0] -
g6[0]/k) == (g12'[0] - g12[0]/k)};


The rest of code is the same and plot for U looks like this one

Update 2. In a case of temperature effect on the flow we have

muw = 0.001004;
muk = 0.00151;
k = 1.5;
\[Rho]Mw = 2100;
\[Rho]Mi = 5060;
\[Rho]Ag = 10500;
\[Rho]w = 997.1;
ph1 = 0.01;
ph2 = 0.01;
ph3 = 0.01;
\[Sigma]Ag = 6.30*10^7;
\[Sigma]w = 5.5*10^-6;
\[Sigma]k = 5*10^-11;
\[Sigma]Mw = 10^-7; \[Sigma]Mi = 2.09*10^4;
btMw = 2.8*10^-5;
btMi = 2.8424*10^-5;
btw = 21*((10)^(-5));
btAg = 5.4*10^-5;
\[Kappa]Ag = 429;
\[Kappa]Mi = 34.5;
\[Kappa]Mw = 3000;
m = 3;
\[Kappa]w = 0.613;
CPMw = 711;
CPw = 4179;
CPMi = 397.746;
CPAg = 235;
CPk = 2090;
\[Rho]k = 783;
\[Kappa]k = 0.15;
btk = 21*((10)^(-5)); M = 1;
b = 5;
l = 2;
Gr = 0.5;
H = 1;
l1 = 2;
rh = \[Rho]w/\[Rho]k;
s = \[Sigma]w/\[Sigma]k;
bt = btw/btk;
M1 = Sqrt[mu/s]*M;
b1 = mu/s*b;
mu = muw/muk;
Gr1 = mu/(rh*bt)*Gr;
H1 = Sqrt[mu/rh]*H;
la = 1;
Pr = 6.8445;
a1 = 1;
a2 = 2;
\[Kappa] = \[Kappa]w/\[Kappa]k;
Pr1 = 21;
x = 0.5;
ep = 0.1;
\[Beta] = 0.001;
t = Pi/3;
A1 = 1/((1 - ph1)^2.5*(1 - ph2)^2.5*(1 - ph3)^2.5);
A2 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw)/\[Rho]w))*(1 - ph2)) + (
ph2*\[Rho]Mi)/\[Rho]w)) + (ph3*\[Rho]Ag)/\[Rho]w;

A3 = ((\[Sigma]Ag*(1 + (2*ph3))) + (2*\[Sigma]w*B3*
C3*(1 - ph3)))/((\[Sigma]Ag*(1 - ph3)) + (\[Sigma]w*C3*
B3*(2 + ph3)));
A4 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*btMw)/(\[Rho]w*btw)))*(1 -
ph2)) + (ph2*\[Rho]Mi*btMi)/(\[Rho]w*btw))) + (
ph3*\[Rho]Ag*btAg)/(\[Rho]w*btw);
A5 = (\[Kappa]Ag + ((m - 1)*C5*B5*\[Kappa]w) - ((m - 1)*
ph3*(C5*B5*\[Kappa]w - \[Kappa]Ag)))/(\[Kappa]Ag + ((m - 1)*C5*
B5*\[Kappa]w) + (ph3*(C5*B5*\[Kappa]w - \[Kappa]Ag)));
A6 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*CPMw)/(\[Rho]w*CPw)))*(1 -
ph2)) + (ph2*\[Rho]Mi*CPMi)/(\[Rho]w*CPw))) + (
ph3*\[Rho]Ag*CPAg)/(\[Rho]w*CPw);
B1 = 1/((1 - ph1)^2.5*(1 - ph2)^2.5);
B2 = (((1 - ph1) + ((ph1*\[Rho]Mw)/\[Rho]k))*(1 - ph2)) + (
ph2*\[Rho]Mi)/\[Rho]k;
B3 = ((\[Sigma]Mi*(1 + (2*ph2))) + (2*\[Sigma]w*
C3*(1 - ph2)))/((\[Sigma]Mi*(1 - ph2)) + (\[Sigma]w*
C3*(2 + ph2)));
B4 = (((1 - ph1) + ((ph1*\[Rho]Mw*btMw)/(btk*\[Rho]k)))*(1 - ph2)) + (
ph2*btMi*\[Rho]Mi)/(btk*\[Rho]k);
B5 = (\[Kappa]Mi + ((m - 1)*C5*\[Kappa]w) - ((m - 1)*
ph2*((C5*\[Kappa]w) - \[Kappa]Mi)))/(\[Kappa]Mi + ((m - 1)*
C5*\[Kappa]w) + (ph2*((C5*\[Kappa]w) - \[Kappa]Mi)));
B6 = (((1 - ph1) + ((ph1*\[Rho]Mw*CPMw)/(\[Rho]k*CPk)))*(1 - ph2)) + (
ph2*\[Rho]Mi*CPMi)/(\[Rho]k*CPk);
C3 = ((\[Sigma]Mw*(1 + (2*ph1))) + (2*\[Sigma]w*(1 -
ph1)))/((\[Sigma]Mw*(1 - ph1)) + (\[Sigma]w*(2 + ph1)));
C5 = (\[Kappa]Mw + ((m - 1)*\[Kappa]w) - ((m - 1)*
ph1*(\[Kappa]w - \[Kappa]Mw)))/(\[Kappa]Mw + ((m -
1)*\[Kappa]w) + (ph1*(\[Kappa]w - \[Kappa]Mw)));
eq = {A1*(u100''[y] + (u100'[y]/(y + k))) - ((A3*B3*C3)*(k/(y + k))^2*
M^2)*u100[y] + ((A3*B3*C3*M*b) + l) + (A4*Gr*T100[y]) == 0,
A1*(u110''[
y] + (u110'[y]/(y + k))) - (((A3*B3*C3)*(k/(y + k))^2*
M^2) + (I*A2*H^2))*u110[y] + l1 + (A4*Gr*T110[y]) == 0,
B1*(u200''[y] + (u200'[y]/(y + k))) - ((B3*C3*s)*(k/(y + k))^2*
M1^2)*u200[y] + ((s*B3*C3*M1*b1) + (mu*l)) + (B4*Gr1*
T200[y]) == 0,
B1*(u210''[
y] + (u210'[y]/(y + k))) - (((B3*C3*s)*(k/(y + k))^2*
M1^2) + (I*B2*H1^2))*u210[y] + (mu*l1) + (B4*Gr1*T210[y]) ==
0, A1*(g1''[
y] + (g1'[y]/(y + k)) - (la^2*(k/(y + k))^2*g1[y])) - ((A3*
B3*C3)*(k/(y + k))^2*M^2)*g1[y] + (A4*Gr*f1[y]) == 0,
A1*(g2''[y] + (g2'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g2[y])) - (((A3*B3*C3)*(k/(y + k))^2*M^2) + (I*A2*H^2))*
g2[y] + (A4*Gr*f2[y]) == 0,
B1*(g3''[y] + (g3'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g3[y])) - ((B3*C3*s)*(k/(y + k))^2*M1^2)*
g3[y] + (B4*Gr1*f3[y]) == 0,
B1*(g4''[
y] + (g4'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g4[y])) - (((B3*C3*s)*(k/(y + k))^2*M1^2) + (I*B2*H1^2))*
g4[y] + (B4*Gr1*f4[y]) == 0,
A1*(m1''[y] + (m1'[y]/(y + k))) - ((A3*B3*C3)*(k/(y + k))^2*M^2)*
m1[y] + (A4*Gr*h1[y]) == 0,
A1*(n1''[
y] + (n1'[y]/(y + k))) - (((A3*B3*C3)*(k/(y + k))^2*
M^2) + (A1*4*la^2*(k/(y + k))^2))*n1[y] + (A4*Gr*k1[y]) ==
0, A1*(m2''[
y] + (m2'[y]/(y + k))) - (((A3*B3*C3)*(k/(y + k))^2*
M^2) + (I*A2*H^2))*m2[y] + (A4*Gr*h2[y]) == 0,
A1*(n2''[
y] + (n2'[
y]/(y + k))) - ((((A3*B3*C3)*(k/(y + k))^2*M^2) + (I*A2*
H^2) + (A1*4*la^2*(k/(y + k))^2)) + (A1*4*
la^2*(k/(y + k))^2))*n2[y] + (A4*Gr*k2[y]) == 0,
B1*(m3''[y] + (m3'[y]/(y + k))) - ((s*B3*C3)*(k/(y + k))^2*M1^2)*
m3[y] + (B4*Gr1*h3[y]) == 0,
B1*(n3''[y] + (n3'[y]/(y + k))) - (((s*B3*C3)*(k/(y + k))^2*
M1^2) + (B1*4*la^2*(k/(y + k))^2))*n3[y] + (B4*Gr1*k3[y]) ==
0,
B1*(m4''[
y] + (m4'[y]/(y + k))) - (((s*B3*C3)*(k/(y + k))^2*
M1^2) + (I*B2*H1^2))*m4[y] + (B4*Gr1*h4[y]) == 0,
B1*(n4''[y] + (n4'[
y]/(y + k))) - ((((s*B3*C3)*(k/(y + k))^2*M1^2) + (I*B2*
H1^2)) + (B1*4*la^2*(k/(y + k))^2))*
n4[y] + (B4*Gr1*h4[y]) == 0,
(((A5*B5*C5)/(A6*Pr))*(T100''[y] + (T100'[y]/(y + k)))) + (a1/A6*
T100[y]) ==
0, (((A5*B5*C5)/(
A6*Pr))*(T110''[y] + (T110'[y]/(y + k)))) + ((a1/
A6 - (I*H^2))*T110[y]) == 0,
(((\[Kappa]*B5*C5)/(
B6*Pr1))*(T200''[y] + (T200'[y]/(y + k)))) + (a2/B6*T200[y]) ==
0, (((\[Kappa]*B5*C5)/(
B6*Pr1))*(T210''[y] + (T200'[y]/(y + k)))) + ((a2/
B6 - (I*H1^2))*T210[y]) ==
0, (((A5*B5*C5)/(
A6*Pr))*(f1''[
y] + (f1'[y]/(y + k)) - (la^2*(k/(y + k))^2*f1[y]))) + (a1/
A6*f1[y]) ==
0, (((A5*B5*C5)/(
A6*Pr))*(f2''[
y] + (f2'[y]/(y + k)) - (la^2*(k/(y + k))^2*f2[y]))) + ((a1/
A6 - (I*H^2))*f2[y]) == 0,
(((\[Kappa]*B5*C5)/(
B6*Pr1))*(f3''[
y] + (f3'[y]/(y + k)) - (la^2*(k/(y + k))^2*f3[y]))) + (a2/
B6*f3[y]) ==
0, (((\[Kappa]*B5*C5)/(
B6*Pr1))*(f4''[
y] + (f4'[y]/(y + k)) - (la^2*(k/(y + k))^2*f4[y]))) + ((a2/
B6 - (I*H1^2))*f4[y]) ==
0, (((A5*B5*C5)/(A6*Pr))*(h1''[y] + (h1'[y]/(y + k)))) + (a1/A6*
h1[y]) == 0,
(((A5*B5*C5)/(
A6*Pr))*(k1''[
y] + (k1'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k1[y]))) + (a1/A6*k1[y]) ==
0, (((A5*B5*C5)/(
A6*Pr))*(h2''[y] + (h2'[y]/(y + k)))) + ((a1/A6 - (I*H^2))*
h2[y]) ==
0, (((A5*B5*C5)/(
A6*Pr))*(k2''[
y] + (k2'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k2[y]))) + ((a1/A6 - (I*H^2))*k2[y]) ==
0, (((\[Kappa]*B5*C5)/(
B6*Pr1))*(h3''[y] + (h3'[y]/(y + k)))) + (a2/B6*h3[y]) ==
0, (((\[Kappa]*B5*C5)/(
B6*Pr1))*(k3''[
y] + (k3'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k3[y]))) + (a2/B6*k3[y]) == 0,
(((\[Kappa]*B5*C5)/(
B6*Pr1))*(h4''[y] + (h4'[y]/(y + k)))) + ((a2/B6 - (I*H1^2))*
h4[y]) ==
0, (((\[Kappa]*B5*C5)/(
B6*Pr1))*(k4''[
y] + (k4'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k4[y]))) + ((a2/B6 - (I*H1^2))*k4[y]) == 0};
bc = {u100[-1] == 0, u110[-1] == 0, u200[1] == 0, u210[1] == 0,
g1[-1] == -u100'[-1], g2[-1] == -u110'[-1], g3[1] == -u200'[1],
g4[1] == -u210'[1], m1[-1] == -(1/2)*(g1'[-1] + (1/2*u100''[-1])),
n1[-1] == (1/2)*(g1'[-1] + (u100''[-1]/2)),
m2[-1] == -(1/2)*(g2'[-1] + (1/2*u110''[-1])),
n2[-1] == (1/2)*(g2'[-1] + (u110''[-1]/2)),
m3[1] == -(1/2)*(g3'[1] + (1/2*u200''[1])),
n3[1] == (1/2)*(g3'[1] + (u200''[1]/2)),
m4[1] == -(1/2)*(g4'[1] + (1/2*u210''[1])),
n4[1] == (1/2)*(g4'[1] + (u210''[1]/2)), u100[0] == u200[0],
u110[0] ==
u210[0], ((mu*A1)/B1)*(u100'[0] - u100[0]/k) == (u200'[0] -
u200[0]/k), ((mu*A1)/B1)*(u110'[0] - u110[0]/k) == (u210'[0] -
u210[0]/k), g1[0] == g3[0],
g2[0] ==
g4[0], ((mu*A1)/B1)*(g1'[0] - g1[0]/k) == (g3'[0] - g3[0]/k), ((
mu*A1)/B1)*(g2'[0] - g2[0]/k) == (g4'[0] - g4[0]/k),
m1[0] == m3[0], m2[0] == m4[0], n1[0] == n3[0],
n2[0] ==
n4[0], ((mu*A1)/B1)*(m1'[0] - m1[0]/k) == (m3'[0] - m3[0]/k), ((
mu*A1)/B1)*(m2'[0] - m2[0]/k) == (m4'[0] - m4[0]/k), ((mu*A1)/
B1)*(n1'[0] - n1[0]/k) == (n3'[0] - n3[0]/k), ((mu*A1)/
B1)*(n2'[0] - n2[0]/k) == (n4'[0] - n4[0]/k), T100[-1] == 0,
T110[-1] == 0, T200[1] == 1, T210[1] == 0, f1[-1] == -T100'[-1],
f2[-1] == -T110'[-1], f3[1] == -T200'[1], f4[1] == -T210'[1],
h1[-1] == (-1/2)*(f1'[-1] + (T100''[-1]/2)),
h2[-1] == (-1/2)*(f2'[-1] + (T110''[-1]/2)),
h3[1] == (-1/2)*(f3'[1] + (T200''[1]/2)),
h4[1] == (-1/2)*(f4'[1] + (T210''[1]/2)),
k1[-1] == (1/2)*(f1'[-1] + (T100''[-1]/2)),
k2[-1] == (1/2)*(f2'[-1] + (T110''[-1]/2)),
k3[1] == (1/2)*(f3'[1] + (T200''[1]/2)),
k4[1] == (1/2)*(f4'[1] + (T210''[1]/2)), T100[0] == T200[0],
f1[0] == f3[0], h1[0] == h3[0], k1[0] == k3[0], T110[0] == T210[0],
f2[0] == f4[0], h2[0] == h4[0], k2[0] == k4[0],
A5*T100'[0] == T200'[0], A5*f1'[0] == f3'[0], A5*h1'[0] == h3'[0],
A5*k1'[0] == k3'[0], A5*T110'[0] == T210'[0], A5*f2'[0] == f4'[0],
A5*h2'[0] == h4'[0], A5*k2'[0] == k4'[0]};
var = {u100, u110, u200, u210, g1, g2, g3, g4, m1, m2, m3, m4, n1, n2,
n3, n4, T100, T110, T200, T210, f1, f2, f3, f4, h1, h2, h3, h4,
k1, k2, k3, k4};
coord = {y, -1, 1};
op = Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}};

var0 = Table[var[[i]]@y, {i, Length[var]}]; var1 = D[var0, y]; var2 =
D[var1, y];

sol0 = Solve[eq, var2];

sol01 = sol0 /. y -> 1; sol02 = sol0 /. y -> -1;

bc1 = bc /. sol01[[1]] /. sol02[[1]];


Numerical solution and visualization

sol = NDSolve[{eq, bc1}, var, coord, op];

Plot[Evaluate[Re[var0] /. sol[[1]]], coord, PlotRange -> All,
PlotLegends -> var, Frame -> True, FrameLabel -> Automatic]


• Respected @AlexTrounev sir, I want my velocity should be zero at 1 and -1 Jun 3, 2023 at 9:10
• @KomalGoyal It sounds good, but what is velocity in your model? Jun 3, 2023 at 11:09
• sir, u10(y,z)+betae^(it)*u11(y,z) is the velocity in one region[0,y,1] sir. And, u20(y,z)+betae^(it)*u21(y,z) is the velocity in another region[-1,y,0] sir Jun 3, 2023 at 12:51
• @KomalGoyal If we put U[-1]==0,U[1]==0 then we need to remove 2 boundary conditions from bc. Could you suggest what we can remove? Jun 3, 2023 at 16:54
• Okay, @Alex Trounev sir, Can you please explain sir, why it is so? Because If we remove any of the boundary conditions from the above equation, then I think we are left with some of the unknown constants. Correct me, please if I am wrong. Jun 4, 2023 at 3:43