1
$\begingroup$

In the present question, I am trying to solve a system of ODEs with corresponding boundary conditions. After that, I tried to find Q, which is firstly dependent on z and then on M. Then, after using numerical integration, I made it independent on z, and now it is dependent on M. Now I am trying to plot Q versus M but I am getting some errors. It might be I made some mistake. Can anyone please check it and help me in that. Your help is highly appreciable Thank you in advance.

Constant numerical values

\[Beta] = 0.1;
    \[Epsilon] = 0.1;
    \[Alpha] = 1;
    \[Zeta] = 0;
    H = I;
    k = 1.5;
    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;

Some variable values we have to find for computations

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]/\[Rho]]*H;

System of Equations

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

Corresponding boundary conditions

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)};

Variables

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

Technique opted to solve the ODEs with boundary conditions

 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]];
    sol = NDSolve[{eq, bc1}, var, coord, op];

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

Some constants using the predefined variables

C1 = Re[NIntegrate[u200[y] /. sol[[1]], {y, -1, 0}]];
C2 = Re[NIntegrate[u100[y] /. sol[[1]], {y, 0, 1}]];
C3 = Re[NIntegrate[u210[y] /. sol[[1]], {y, -1, 0}]];
C4 = Re[NIntegrate[u110[y] /. sol[[1]], {y, 0, 1}]];
C5 = Re[u200[y] /. sol[[1]] /. y -> 0];
C6 = Re[u100[y] /. sol[[1]] /. y -> 0];
C7 = Re[u210[y] /. sol[[1]] /. y -> 0];
C8 = Re[u110[y] /. sol[[1]] /. y -> 0];
C9[(z_)?NumericQ] := NIntegrate[u201[y, z], {y, -1, 0}];
C10[(z_)?NumericQ] := NIntegrate[u101[y, z], {y, 0, 1}];
C11[(z_)?NumericQ] := NIntegrate[u211[y, z], {y, -1, 0}];
C12[(z_)?NumericQ] := NIntegrate[u111[y, z], {y, 0, 1}];
C13 = Re[D[u200[y], y] /. sol[[1]] /. y -> 0];
C14 = Re[D[u100[y], y] /. sol[[1]] /. y -> 0];
C15 = Re[D[u110[y], y] /. sol[[1]] /. y -> 0];
C16 = Re[D[u210[y], y] /. sol[[1]] /. y -> 0];
C17[z_?NumericQ] := u201[y, z] /. y -> 0;
C18[z_?NumericQ] := u201[y, z] /. y -> -1;
C19[z_?NumericQ] := u101[y, z] /. y -> 1;
C20[z_?NumericQ] := u101[y, z] /. y -> 0;
C21[z_?NumericQ] := u211[y, z] /. y -> 0;
C22[z_?NumericQ] := u211[y, z] /. y -> -1;
C23[z_?NumericQ] := u111[y, z] /. y -> 1;
C24[z_?NumericQ] := u111[y, z] /. y -> 0;
C25[z_?NumericQ] := NIntegrate[u202[y, z], {y, -1, 0}];
C26[z_?NumericQ] := NIntegrate[u102[y, z], {y, 0, 1}];
C27[z_?NumericQ] := NIntegrate[u212[y, z], {y, -1, 0}];
C28[z_?NumericQ] := NIntegrate[u112[y, z], {y, 0, 1}];

Numerical intergral

R[z_, M_] := 
 Re[(C1 + C2) + (\[Beta]*
     E^(I*t)*(C3 + 
       C4)) + (\[Epsilon]*((-Sin[\[Alpha]*z + \[Zeta]]*(C5 + (\[Beta]*
             E^(I*t)*C7))) + (-Sin[\[Alpha]*
            z]*(C6 + (\[Beta]*E^(I*t)*C8))) + (C9[z] + 
         C10[z]) + ((\[Beta]*E^(I*t))*(C11[z] + 
           C12[z])))) + (\[Epsilon]^2*((-C13*
          Sin[\[Alpha]*z + \[Zeta]]^2/2 - 
         C14*Sin[\[Alpha]*z]^2/2) + (\[Beta]*
         E^(I*t)*(-C16*Sin[\[Alpha]*z + \[Zeta]]^2/2 - 
           C15*Sin[\[Alpha]*z]^2/
             2)) + ((-Sin[\[Alpha]*z + \[Zeta]]*(C17[z] - 
             C18[z])) + (Sin[\[Alpha]*z]*(C19[z] - 
             C20[z]))) + (\[Beta]*
         E^(I*t)*((-Sin[\[Alpha]*z + \[Zeta]]*(C21[z] - 
               C22[z])) + (Sin[\[Alpha]*z]*(C23[z] - 
               C24[z])))) + (C25[z] + C26[z]) + (\[Beta]*
         E^(I*t)*(C27[z] + C28[z]))))]

Q dependent on M

Q[M_] := (\[Alpha] Re[
         NIntegrate[
          Evaluate[R[z, M]], {z, 0, (2 \[Pi])/\[Alpha]}]])/(2 \[Pi])

Final Plot

P1 = Plot[Q[M], {M, 0, 10}, 
  PlotLegends -> 
   Placed[LineLegend[{Style[
       "\!\(\*SubscriptBox[\(M\), \(Water\)]\)=0", Bold, 12]}], {0.2, 
     0.2}], PlotStyle -> {Cyan}, PlotRange -> All, Axes -> False, 
  Frame -> True, 
  FrameLabel -> {Style["t", Black, 13, Bold], 
    Style["Q", Black, 13, Bold]}, 
  LabelStyle -> 
   Directive[{Black, 15}, Bold, FontFamily -> "Times new Roman"]]
$\endgroup$
8
  • $\begingroup$ @AlexTrounev, sir, can you please help me in finding that plot against M. I tried but I think I am doing some mistake. $\endgroup$ Commented Jul 2, 2023 at 5:03
  • $\begingroup$ As I remember we solve this in previous post but then post with answer has been deleted. Also my answer has not been accepted on mathematica.stackexchange.com/questions/286292/… :) $\endgroup$ Commented Jul 2, 2023 at 5:56
  • $\begingroup$ @AlexTrounev Yes, sir, because some mistake was there in that code, I want to give the reference I was uploading that code which create a problem for me. Hence I was deleting some more code also even together with that one. And, at the time you answered my question, I don't know that I had to accept answers. Now it's done sir I accept that answer. You can check now. $\endgroup$ Commented Jul 2, 2023 at 6:17
  • $\begingroup$ Thank you. There is a typo in your code. Parameter M is not defined. $\endgroup$ Commented Jul 2, 2023 at 7:51
  • $\begingroup$ @AlexTrounev how to define M sir, If I am giving some value to M then it is showing so many errors since I want M on x-axis of the plot and Q on the y-axis $\endgroup$ Commented Jul 2, 2023 at 9:26

1 Answer 1

1
$\begingroup$

At M>7 solution is not stable, therefore we plot Q[M] in the range $0\le M\le 7$ as follows

\[Beta] = 0.1;
\[Epsilon] = 0.1;
\[Alpha] = 1;
\[Zeta] = 0;
H = I;
k = 1.5;
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;
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]/\[Rho]]*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]];
sol1[MM_] := NDSolve[{eq, bc1} /. {M -> MM}, var, coord, op];

Do[sol = sol1[j][[1]];
   With[{n1 = sol}, 
    u101[y_, z_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. 
      n1;
    u111[y_, z_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. 
      n1;
    u201[y_, z_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. 
      n1;
    u211[y_, z_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. 
      n1;
    u102[y_, 
      z_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + 
        g3[y]*Cos[2*\[Alpha]*z]) /. n1;
    u112[y_, 
      z_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + 
        g6[y]*Cos[2*\[Alpha]*z]) /. n1;
    u202[y_, 
      z_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + 
        g9[y]*Cos[2*\[Alpha]*z]) /. n1;
    u212[y_, 
      z_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + 
        g12[y]*Cos[2*\[Alpha]*z]) /. n1;
    u10[y_, 
      z_] := (u100[y] + \[Epsilon]*u101[y, z] + \[Epsilon]^2*
         u102[y, z]) /. n1;
    u11[y_, 
      z_] := (u110[y] + \[Epsilon]*u111[y, z] + \[Epsilon]^2*
         u112[y, z]) /. n1;
    u20[y_, 
      z_] := (u200[y] + \[Epsilon]*u201[y, z] + \[Epsilon]^2*
         u202[y, z]) /. n1;
    u21[y_, 
      z_] := (u210[y] + \[Epsilon]*u211[y, z] + \[Epsilon]^2*
         u212[y, z]) /. n1]; 
   C17[z_?NumericQ] := u201[y, z] /. y -> 0;
   C18[z_?NumericQ] := u201[y, z] /. y -> -1;
   C19[z_?NumericQ] := u101[y, z] /. y -> 1;
   C20[z_?NumericQ] := u101[y, z] /. y -> 0;
   C21[z_?NumericQ] := u211[y, z] /. y -> 0;
   C22[z_?NumericQ] := u211[y, z] /. y -> -1;
   C23[z_?NumericQ] := u111[y, z] /. y -> 1;
   C24[z_?NumericQ] := u111[y, z] /. y -> 0;
   C25[z_?NumericQ] := NIntegrate[u202[y, z], {y, -1, 0}];
   C26[z_?NumericQ] := NIntegrate[u102[y, z], {y, 0, 1}];
   C27[z_?NumericQ] := NIntegrate[u212[y, z], {y, -1, 0}];
   C28[z_?NumericQ] := NIntegrate[u112[y, z], {y, 0, 1}];
   
   C9[(z_)?NumericQ] := NIntegrate[u201[y, z], {y, -1, 0}];
   C10[(z_)?NumericQ] := NIntegrate[u101[y, z], {y, 0, 1}];
   C11[(z_)?NumericQ] := NIntegrate[u211[y, z], {y, -1, 0}];
   C12[(z_)?NumericQ] := NIntegrate[u111[y, z], {y, 0, 1}];
   
   C1 = Re[NIntegrate[u200[y] /. sol, {y, -1, 0}]];
   C2 = Re[NIntegrate[u100[y] /. sol, {y, 0, 1}]];
   C3 = Re[NIntegrate[u210[y] /. sol, {y, -1, 0}]];
   C4 = Re[NIntegrate[u110[y] /. sol, {y, 0, 1}]];
   C5 = Re[u200[y] /. sol /. y -> 0];
   C6 = Re[u100[y] /. sol /. y -> 0];
   C7 = Re[u210[y] /. sol /. y -> 0];
   C8 = Re[u110[y] /. sol /. y -> 0];
   
   C13 = Re[D[u200[y], y] /. sol /. y -> 0];
   C14 = Re[D[u100[y], y] /. sol /. y -> 0];
   C15 = Re[D[u110[y], y] /. sol /. y -> 0];
   C16 = Re[D[u210[y], y] /. sol /. y -> 0]; 
   R[t_, z_, M_] := 
    Re[(C1 + C2) + (\[Beta]*
        E^(I*t)*(C3 + 
          C4)) + (\[Epsilon]*((-Sin[\[Alpha]*
                z + \[Zeta]]*(C5 + (\[Beta]*E^(I*t)*
                C7))) + (-Sin[\[Alpha]*
               z]*(C6 + (\[Beta]*E^(I*t)*C8))) + (C9[z] + 
            C10[z]) + ((\[Beta]*E^(I*t))*(C11[z] + 
              C12[z])))) + (\[Epsilon]^2*((-C13*
             Sin[\[Alpha]*z + \[Zeta]]^2/2 - 
            C14*Sin[\[Alpha]*z]^2/2) + (\[Beta]*
            E^(I*t)*(-C16*Sin[\[Alpha]*z + \[Zeta]]^2/2 - 
              C15*Sin[\[Alpha]*z]^2/
                2)) + ((-Sin[\[Alpha]*z + \[Zeta]]*(C17[z] - 
                C18[z])) + (Sin[\[Alpha]*z]*(C19[z] - 
                C20[z]))) + (\[Beta]*
            E^(I*t)*((-Sin[\[Alpha]*z + \[Zeta]]*(C21[z] - 
                  C22[z])) + (Sin[\[Alpha]*z]*(C23[z] - 
                  C24[z])))) + (C25[z] + C26[z]) + (\[Beta]*
            E^(I*t)*(C27[z] + C28[z]))))]; 
   Q[j] = (\[Alpha] Re[
        NIntegrate[
         Evaluate[R[1, z, j]], {z, 
          0, (2 \[Pi])/\[Alpha]}]])/(2 \[Pi]);, {j, 0, 7, 
    1/2}];

It takes a time, but finally we have

QQ = Interpolation[Table[{j, Q[j]}, {j, 0, 7, 1/2}]];

Plot[QQ[M], {M, 0, 7}, FrameLabel -> {"M", "Q"}, Frame -> True]

Figure 1

Note, that we plot Q[M] at t=1, if you need data for t=Pi/3, then use in the last line Do loop

Q[j] = (\[Alpha] Re[
            NIntegrate[
             Evaluate[R[Pi/3, z, j]], {z, 
              0, (2 \[Pi])/\[Alpha]}]])/(2 \[Pi]);
$\endgroup$
2
  • $\begingroup$ Thank you @AlexTrounev, for your help can you suggest if we want to fix M and want k on the x-axis. How do we proceed because k is in the boundary condition? So I am facing a little difficulty while plotting Q against k $\endgroup$ Commented Jul 3, 2023 at 14:19
  • $\begingroup$ Do you mean Q[k] at fixed M? First define M, for instance, M=1;, then define sol1[MM_] := NDSolve[{eq, bc1} /. {k -> MM}, var, coord, op];. Finally use code, compute Q[j] and plot QQ[k]. $\endgroup$ Commented Jul 4, 2023 at 2:20

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.