1
$\begingroup$

enter image description hereIn 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.

Thanking You in Advance

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}]
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Kuba
    Jun 13, 2023 at 6:26

1 Answer 1

3
$\begingroup$

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]

Figure 1

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}]

Figure 2

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 Figure 3

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]

Figure 4

$\endgroup$
23
  • $\begingroup$ Respected @AlexTrounev sir, I want my velocity should be zero at 1 and -1 $\endgroup$ Jun 3, 2023 at 9:10
  • $\begingroup$ @KomalGoyal It sounds good, but what is velocity in your model? $\endgroup$ Jun 3, 2023 at 11:09
  • $\begingroup$ 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 $\endgroup$ Jun 3, 2023 at 12:51
  • $\begingroup$ @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? $\endgroup$ Jun 3, 2023 at 16:54
  • $\begingroup$ 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. $\endgroup$ Jun 4, 2023 at 3:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.