# Plotting of graph after solving system of ODEs using NDSolve command

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

• @AlexTrounev, sir, can you please help me in finding that plot against M. I tried but I think I am doing some mistake. Commented Jul 2, 2023 at 5:03
• 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/… :) Commented Jul 2, 2023 at 5:56
• @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. Commented Jul 2, 2023 at 6:17
• Thank you. There is a typo in your code. Parameter M is not defined. Commented Jul 2, 2023 at 7:51
• @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 Commented Jul 2, 2023 at 9:26

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]


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

• 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 Commented Jul 3, 2023 at 14:19
• 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]. Commented Jul 4, 2023 at 2:20