# R-K fourth order scheme for system of ODEs containing 32 differential equations

In the present code, I am trying to solve the system of differential equations using the R-K fourth-order scheme. I am trying to plot the function u and T, but it is showing some error in the boundary conditions. Can anyone please help to identify the error? Suppose you need any further information. Kindly drop a comment. Thank you in advance.

Values for the constants

ph1 = 0.01; ph2 = 0.01; ph3 = 0.01; \[Rho]Mw = 2100; \[Rho]Mi = 5060; \
\[Rho]Ag = 10500; \[Rho]w = 997.1;
\[Rho]k = 783; \[Sigma]Mw = 10^-7; \[Sigma]Mi =
2.09*10^4; \[Sigma]Ag = 6.30*10^7; \[Sigma]w = 5.5*10^-6; \[Sigma]k =
5*10^-11;
btMw = 2.8*10^-5; btMi = 2.8424*10^-5; btAg = 5.4*10^-5; btw =
21*((10)^(-5));
btk = 21*((10)^(-5)); CPMw = 711; CPMi = 397.746; CPAg = 235; CPw = \
4179; CPk = 2090;
\[Kappa]Mw = 3000; \[Kappa]Mi = 34.5; \[Kappa]Ag = 429; \[Kappa]w = \
0.613; \[Kappa]k = 0.15; muw = 0.001004; muk = 0.00151;


Ratios for the constants

mu = muw/muk; \[Rho] = \[Rho]w/\[Rho]k; \[Sigma] = \
\[Sigma]w/\[Sigma]k; bt = btw/btk; CP = CPw/CPk; \[Kappa] = \
\[Kappa]w/\[Kappa]k;


Physical parameters

m = 3; l = 3; l1 = 2; ep = 0.1; \[Beta] = 0.001; t =
Pi/3; \[Beta]1 = 1; \[Alpha] = 0; k = 1.5; la = 1; x = 10;


Some parameters

Gr = (Gr1*\[Rho]*bt)/mu; M = M1*Sqrt[\[Sigma]/mu]; b =
b1*Sqrt[\[Sigma]/mu]; H =
H1*Sqrt[\[Rho]/mu]; Pr = 6.84; \[CapitalOmega] = \[CapitalOmega]1/(
mu*CP);


Some more constants

M1 = 1; b1 = 5; Gr1 = 0.5; H1 = 1; Pr1 = 21; \[CapitalOmega]1 = 1;


Some more ratios

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)));
B3 = ((\[Sigma]Mi*(1 + (2*ph2))) + (2*\[Sigma]w*
C3*(1 - ph2)))/((\[Sigma]Mi*(1 - ph2)) + (\[Sigma]w*
C3*(2 + ph2)));
C3 = ((\[Sigma]Mw*(1 + (2*ph1))) + (2*\[Sigma]w*(1 -
ph1)))/((\[Sigma]Mw*(1 - ph1)) + (\[Sigma]w*(2 + ph1)));
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 = ((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);
A6 = (\[Kappa]Ag + ((m - 1)*C6*B6*\[Kappa]w) - ((m - 1)*
ph3*(C6*B6*\[Kappa]w - \[Kappa]Ag)))/(\[Kappa]Ag + ((m - 1)*C6*
B6*\[Kappa]w) + (ph3*(C6*B6*\[Kappa]w - \[Kappa]Ag)));
B6 = (\[Kappa]Mi + ((m - 1)*C6*\[Kappa]w) - ((m - 1)*
ph2*((C6*\[Kappa]w) - \[Kappa]Mi)))/(\[Kappa]Mi + ((m - 1)*
C6*\[Kappa]w) + (ph2*((C6*\[Kappa]w) - \[Kappa]Mi)));
C6 = (\[Kappa]Mw + ((m - 1)*\[Kappa]w) - ((m - 1)*
ph1*(\[Kappa]w - \[Kappa]Mw)))/(\[Kappa]Mw + ((m -
1)*\[Kappa]w) + (ph1*(\[Kappa]w - \[Kappa]Mw)));

A11 = 1/((1 - ph1)^2.5*(1 - ph2)^2.5*(1 - ph3)^2.5);
A12 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw)/\[Rho]k))*(1 - ph2)) + (
ph2*\[Rho]Mi)/\[Rho]k)) + (ph3*\[Rho]Ag)/\[Rho]k;
A13 = ((\[Sigma]Ag*(1 + (2*ph3))) + (2*\[Sigma]k*B13*
C13*(1 - ph3)))/((\[Sigma]Ag*(1 - ph3)) + (\[Sigma]k*C13*
B13*(2 + ph3)));
B13 = ((\[Sigma]Mi*(1 + (2*ph2))) + (2*\[Sigma]k*
C13*(1 - ph2)))/((\[Sigma]Mi*(1 - ph2)) + (\[Sigma]k*
C13*(2 + ph2)));
C13 = ((\[Sigma]Mw*(1 + (2*ph1))) + (2*\[Sigma]k*(1 -
ph1)))/((\[Sigma]Mw*(1 - ph1)) + (\[Sigma]k*(2 + ph1)));
A14 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*btMw)/(\[Rho]k*btk)))*(1 -
ph2)) + (ph2*\[Rho]Mi*btMi)/(\[Rho]k*btk))) + (
ph3*\[Rho]Ag*btAg)/(\[Rho]k*btk);
A15 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*CPMw)/(\[Rho]k*CPk)))*(1 -
ph2)) + (ph2*\[Rho]Mi*CPMi)/(\[Rho]k*CPk))) + (
ph3*\[Rho]Ag*CPAg)/(\[Rho]k*CPk);
A16 = (\[Kappa]Ag + ((m - 1)*C16*B16*\[Kappa]k) - ((m - 1)*
ph3*(C16*B16*\[Kappa]k - \[Kappa]Ag)))/(\[Kappa]Ag + ((m - 1)*
C16*B16*\[Kappa]k) + (ph3*(C16*B16*\[Kappa]k - \[Kappa]Ag)));
B16 = (\[Kappa]Mi + ((m - 1)*C16*\[Kappa]k) - ((m - 1)*
ph2*((C16*\[Kappa]k) - \[Kappa]Mi)))/(\[Kappa]Mi + ((m - 1)*
C16*\[Kappa]k) + (ph2*((C16*\[Kappa]k) - \[Kappa]Mi)));
C16 = (\[Kappa]Mw + ((m - 1)*\[Kappa]k) - ((m - 1)*
ph1*(\[Kappa]k - \[Kappa]Mw)))/(\[Kappa]Mw + ((m -
1)*\[Kappa]k) + (ph1*(\[Kappa]k - \[Kappa]Mw)));


Equations

With[{soln = {1 -> 1}}, u101[y_] := Sin[la*x]*g1[y] /. soln;
u102[y_] := (m1[y] + Cos[2*la*x]*n1[y]) /. soln;
u111[y_] := Sin[la*x]*g2[y] /. soln;
u112[y_] := (m2[y] + Cos[2*la*x]*n2[y]) /. soln;
u201[y_] := Sin[la*x]*g3[y] /. soln;
u202[y_] := (m3[y] + Cos[2*la*x]*n3[y]) /. soln;
u211[y_] := Sin[la*x]*g4[y] /. soln;
u212[y_] := (m4[y] + Cos[2*la*x]*n4[y]) /. soln;
u10[y_] := (u100[y] + ep*u101[y] + ep^2*u102[y]) /. soln;
u11[y_] := (u110[y] + ep*u111[y] + ep^2*u112[y]) /. soln;
u20[y_] := (u200[y] + ep*u201[y] + ep^2*u202[y]) /. soln;
u21[y_] := (u210[y] + ep*u211[y] + ep^2*u212[y]) /. soln;
T101[y_] := Sin[la*x]*f1[y] /. soln;
T102[y_] := (h1[y] + Cos[2*la*x]*k1[y]) /. soln;
T111[y_] := Sin[la*x]*f2[y] /. soln;
T112[y_] := (h2[y] + Cos[2*la*x]*k2[y]) /. soln;
T201[y_] := Sin[la*x]*f3[y] /. soln;
T202[y_] := (h3[y] + Cos[2*la*x]*k3[y]) /. soln;
T211[y_] := Sin[la*x]*f4[y] /. soln;
T212[y_] := (h4[y] + Cos[2*la*x]*k4[y]) /. soln;
T10[y_] := (T100[y] + ep*T101[y] + ep^2*T102[y]) /. soln;
T11[y_] := (T110[y] + ep*T111[y] + ep^2*T112[y]) /. soln;
T20[y_] := (T200[y] + ep*T201[y] + ep^2*T202[y]) /. soln;
T21[y_] := (T210[y] + ep*T211[y] + ep^2*T212[y]) /. soln];
U[x_, y_] :=
Piecewise[{{(u20[y] + (\[Beta]*(E^(I*t))*u21[y])), \[Alpha] <= y <=
1}, {(u10[y] + (\[Beta]*(E^(I*t))*u11[y])), -1 <= y <= \[Alpha]}}]
eq = {A1*(1 +
1/\[Beta]1)*(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*(1 + 1/\[Beta]1)*(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,
A11*(u200''[y] + (u200'[y]/(y + k))) - ((A13*B13*C13)*(k/(
y + k))^2*M1^2)*
u200[y] + ((A13*B13*C13*M1*b1) + (mu*l)) + (A14*Gr1*T200[y]) ==
0, A11*(u210''[
y] + (u210'[y]/(y + k))) - (((A13*B13*C13)*(k/(y + k))^2*
M1^2) + (I*A12*H1^2))*
u210[y] + (mu*l1) + (A14*Gr1*T210[y]) == 0,
A1*(1 + 1/\[Beta]1)*(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*(1 + 1/\[Beta]1)*(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,
A11*(g3''[
y] + (g3'[y]/(y + k)) - (la^2*(k/(y + k))^2*g3[y])) - ((A13*
B13*C13)*(k/(y + k))^2*M1^2)*g3[y] + (A14*Gr1*f3[y]) == 0,
A11*(g4''[
y] + (g4'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g4[y])) - (((A13*B13*C13)*(k/(y + k))^2*M1^2) + (I*A12*
H1^2))*g4[y] + (A14*Gr1*f4[y]) == 0,
A1*(1 +
1/\[Beta]1)*(m1''[y] + (m1'[y]/(y + k))) - ((A3*B3*C3)*(k/(
y + k))^2*M^2)*m1[y] + (A4*Gr*h1[y]) == 0,
A1*(1 +
1/\[Beta]1)*(n1''[
y] + (n1'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*n1[y])) - ((A3*
B3*C3)*(k/(y + k))^2*M^2)*n1[y] + (A4*Gr*k1[y]) == 0,
A1*(1 +
1/\[Beta]1)*(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*(1 +
1/\[Beta]1)*(n2''[
y] + (n2'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n2[y])) - (((A3*B3*C3)*(k/(y + k))^2*M^2) + (I*A2*H^2))*
n2[y] + (A4*Gr*k2[y]) == 0,
A11*(m3''[y] + (m3'[y]/(y + k))) - ((A13*B13*C13)*(k/(y + k))^2*
M1^2)*m3[y] + (A14*Gr1*h3[y]) == 0,
A11*(n3''[y] + (n3'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n3[y])) - ((A13*B13*C13)*(k/(y + k))^2*M1^2)*
n3[y] + (A14*Gr1*k3[y]) == 0,
A11*(m4''[
y] + (m4'[y]/(y + k))) - (((A13*B13*C13)*(k/(y + k))^2*
M1^2) + (I*A12*H1^2))*m4[y] + (A14*Gr1*h4[y]) == 0,
A11*(n4''[y] + (n4'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n4[y])) - (((A13*B13*C13)*(k/(y + k))^2*M1^2) + (I*A12*
H1^2))*n4[y] + (A14*Gr1*k4[y]) == 0,
(((A6*B6*C6)/
Pr)*(T100''[y] + (T100'[y]/(y + k)))) + (\[CapitalOmega]*
T100[y]) ==
0, (((A6*B6*C6)/
Pr)*(T110''[
y] + (T110'[y]/(y + k)))) + ((\[CapitalOmega] - (I*A5*H^2))*
T110[y]) == 0,
(((A16*B16*C16)/
Pr1)*(T200''[y] + (T200'[y]/(y + k)))) + (\[CapitalOmega]1*
T200[y]) ==
0, (((A16*B16*C16)/
Pr1)*(T210''[
y] + (T200'[y]/(y + k)))) + ((\[CapitalOmega]1 - (I*A15*
H1^2))*T210[y]) ==
0, (((A6*B6*C6)/
Pr)*(f1''[y] + (f1'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f1[y]))) + (\[CapitalOmega]*f1[y]) ==
0, (((A6*B6*C6)/
Pr)*(f2''[y] + (f2'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f2[y]))) + ((\[CapitalOmega] - (I*A5*H^2))*f2[y]) == 0,
(((A16*B16*C16)/
Pr1)*(f3''[
y] + (f3'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f3[y]))) + (\[CapitalOmega]1*f3[y]) ==
0, (((A16*B16*C16)/
Pr1)*(f4''[
y] + (f4'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f4[y]))) + ((\[CapitalOmega]1 - (I*A15*H1^2))*f4[y]) ==
0, (((A6*B6*C6)/
Pr)*(h1''[y] + (h1'[y]/(y + k)))) + (\[CapitalOmega]*h1[y]) ==
0,
(((A6*B6*C6)/
Pr)*(k1''[y] + (k1'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k1[y]))) + (\[CapitalOmega]*k1[y]) ==
0, (((A6*B6*C6)/
Pr)*(h2''[y] + (h2'[y]/(y + k)))) + ((\[CapitalOmega] - (I*A5*
H^2))*h2[y]) ==
0, (((A6*B6*C6)/
Pr)*(k2''[y] + (k2'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k2[y]))) + ((\[CapitalOmega] - (I*A5*H^2))*k2[y]) ==
0, (((A16*B16*C16)/
Pr1)*(h3''[y] + (h3'[y]/(y + k)))) + (\[CapitalOmega]1*
h3[y]) ==
0, (((A16*B16*C16)/
Pr1)*(k3''[
y] + (k3'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k3[y]))) + (\[CapitalOmega]1*k3[y]) == 0,
(((A16*B16*C16)/
Pr1)*(h4''[
y] + (h4'[y]/(y + k)))) + ((\[CapitalOmega]1 - (I*A15*
H1^2))*h4[y]) ==
0, (((A16*B16*C16)/
Pr1)*(k4''[
y] + (k4'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k4[y]))) + ((\[CapitalOmega]1 - (I*A15*H1^2))*k4[y]) ==
0};


Boundary conditions

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[\[Alpha]] == u200[\[Alpha]],
u110[\[Alpha]] ==
u210[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(u100'[\[Alpha]]) == (u200'[\[Alpha]]), (mu*(1 \
+ 1/\[Beta]1))*(u110'[\[Alpha]]) == (u210'[\[Alpha]]),
g1[\[Alpha]] == g3[\[Alpha]],
g2[\[Alpha]] ==
g4[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(g1'[\[Alpha]]) == (g3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(g2'[\[Alpha]]) == (g4'[\[Alpha]]),
m1[\[Alpha]] == m3[\[Alpha]], m2[\[Alpha]] == m4[\[Alpha]],
n1[\[Alpha]] == n3[\[Alpha]],
n2[\[Alpha]] ==
n4[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(m1'[\[Alpha]]) == (m3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(m2'[\[Alpha]]) == (m4'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(n1'[\[Alpha]]) == (n3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(n2'[\[Alpha]]) == (n4'[\[Alpha]]),
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[\[Alpha]] == T200[\[Alpha]], f1[\[Alpha]] == f3[\[Alpha]],
h1[\[Alpha]] == h3[\[Alpha]], k1[\[Alpha]] == k3[\[Alpha]],
T110[\[Alpha]] == T210[\[Alpha]], f2[\[Alpha]] == f4[\[Alpha]],
h2[\[Alpha]] == h4[\[Alpha]],
k2[\[Alpha]] ==
k4[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
T100'[\[Alpha]] ==
T200'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
f1'[\[Alpha]] ==
f3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*h1'[\[Alpha]] ==
h3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
k1'[\[Alpha]] ==
k3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
T110'[\[Alpha]] ==
T210'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
f2'[\[Alpha]] ==
f4'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*h2'[\[Alpha]] ==
h4'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
k2'[\[Alpha]] == k4'[\[Alpha]]};


Solution Methodology

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]];
sol = NDSolve[{eq, bc1}, var, coord, op];
Plot[Evaluate[Re[var0] /. sol[[1]]], coord, PlotRange -> All,
PlotLegends -> var, Frame -> True, FrameLabel -> Automatic]
{p1 = Plot[U[y] /. sol[[1]] // Re, {y, -1, 1}, PlotRange -> All],
p2 = Plot[
Piecewise[{{Re[(T20[y] + (\[Beta]*(E^(I*t))*T21[y]))], \[Alpha] <=
y <= 1}, {Re[(T10[y] + (\[Beta]*(E^(I*t))*T11[y]))], -1 <=
y <= \[Alpha]}}] /. sol[[1]], {y, -1, 1}, PlotRange -> All]}

• There is messege from NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly. Here you again use the 2nd derivative in the boundary conditions. This is prohibited in the solution methods used by the NDSolve. See this post mathematica.stackexchange.com/questions/285960/… Commented Oct 1, 2023 at 6:26
• @AlexTrounev Okay, sir, I will edit my accordingly. Thank you so much for the clarity sir Commented Oct 1, 2023 at 6:45
• Please, check this parameter ((A3*B3*C3*M*b) + l) in eq[[1]]. Now it is about 865333. Also parameter M^2 is about 165438. It is too high for numerical integration. Commented Oct 1, 2023 at 7:44
• @AlexTrounev sir, if I am not wrong, in mathematica.stackexchange.com/questions/285960/… we try to plot velocity and temperature separately, but in the present case we cannot do that, the equation is coupled in velocity equation temperature is also coming like h1, h2, h3 .... and so onn. And if it taking such a big value we can plot for M=0 also, sir. Commented Oct 1, 2023 at 8:36
• Then see my answer for M^2=1.5. Commented Oct 1, 2023 at 11:29

Parameter M^2 is about 165438. It is too high for numerical integration. If we reduce this parameter to M^2=1.5 then we have

(*Values for the constants*)

ph1 = 0.01; ph2 = 0.01; ph3 = 0.01; \[Rho]Mw = 2100; \[Rho]Mi = 5060; \
\[Rho]Ag = 10500; \[Rho]w = 997.1;
\[Rho]k = 783; \[Sigma]Mw = 10^-7; \[Sigma]Mi =
2.09*10^4; \[Sigma]Ag = 6.30*10^7; \[Sigma]w =
5.5*10^-6; \[Sigma]k = \[Sigma]w(*5*10^-11*);
btMw = 2.8*10^-5; btMi = 2.8424*10^-5; btAg = 5.4*10^-5; btw =
21*((10)^(-5));
btk = 21*((10)^(-5)); CPMw = 711; CPMi = 397.746; CPAg = 235; CPw = \
4179; CPk = 2090;
\[Kappa]Mw = 3000; \[Kappa]Mi = 34.5; \[Kappa]Ag = 429; \[Kappa]w = \
0.613; \[Kappa]k = 0.15; muw = 0.001004; muk = 0.00151;
(*Ratios for the constants*)

mu = muw/
muk; \[Rho] = \[Rho]w/\[Rho]k; \[Sigma] = \[Sigma]w/\[Sigma]k; bt =
btw/btk; CP = CPw/CPk; \[Kappa] = \[Kappa]w/\[Kappa]k;
(*Physical parameters*)

m = 3; l = 3; l1 = 2; ep = 0.1; \[Beta] = 0.001; t =
Pi/3; \[Beta]1 = 1; \[Alpha] = 0; k = 1.5; la = 1; x = 10;
(*Some parameters*)

Gr = (Gr1*\[Rho]*bt)/mu; M = M1*Sqrt[\[Sigma]/mu]; b =
b1*Sqrt[\[Sigma]/mu]; H =
H1*Sqrt[\[Rho]/
mu]; Pr = 6.84; \[CapitalOmega] = \[CapitalOmega]1/(mu*CP);
(*Some more constants*)

M1 = 1; b1 = 5; Gr1 = 0.5; H1 = 1; Pr1 = 21; \[CapitalOmega]1 = 1;
(*Some more ratios*)

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)));
B3 = ((\[Sigma]Mi*(1 + (2*ph2))) + (2*\[Sigma]w*
C3*(1 - ph2)))/((\[Sigma]Mi*(1 - ph2)) + (\[Sigma]w*
C3*(2 + ph2)));
C3 = ((\[Sigma]Mw*(1 + (2*ph1))) + (2*\[Sigma]w*(1 -
ph1)))/((\[Sigma]Mw*(1 - ph1)) + (\[Sigma]w*(2 + ph1)));
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 = ((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);
A6 = (\[Kappa]Ag + ((m - 1)*C6*B6*\[Kappa]w) - ((m - 1)*
ph3*(C6*B6*\[Kappa]w - \[Kappa]Ag)))/(\[Kappa]Ag + ((m - 1)*C6*
B6*\[Kappa]w) + (ph3*(C6*B6*\[Kappa]w - \[Kappa]Ag)));
B6 = (\[Kappa]Mi + ((m - 1)*C6*\[Kappa]w) - ((m - 1)*
ph2*((C6*\[Kappa]w) - \[Kappa]Mi)))/(\[Kappa]Mi + ((m - 1)*
C6*\[Kappa]w) + (ph2*((C6*\[Kappa]w) - \[Kappa]Mi)));
C6 = (\[Kappa]Mw + ((m - 1)*\[Kappa]w) - ((m - 1)*
ph1*(\[Kappa]w - \[Kappa]Mw)))/(\[Kappa]Mw + ((m -
1)*\[Kappa]w) + (ph1*(\[Kappa]w - \[Kappa]Mw)));

A11 = 1/((1 - ph1)^2.5*(1 - ph2)^2.5*(1 - ph3)^2.5);
A12 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw)/\[Rho]k))*(1 -
ph2)) + (ph2*\[Rho]Mi)/\[Rho]k)) + (ph3*\[Rho]Ag)/\[Rho]k;
A13 = ((\[Sigma]Ag*(1 + (2*ph3))) + (2*\[Sigma]k*B13*
C13*(1 - ph3)))/((\[Sigma]Ag*(1 - ph3)) + (\[Sigma]k*C13*
B13*(2 + ph3)));
B13 = ((\[Sigma]Mi*(1 + (2*ph2))) + (2*\[Sigma]k*
C13*(1 - ph2)))/((\[Sigma]Mi*(1 - ph2)) + (\[Sigma]k*
C13*(2 + ph2)));
C13 = ((\[Sigma]Mw*(1 + (2*ph1))) + (2*\[Sigma]k*(1 -
ph1)))/((\[Sigma]Mw*(1 - ph1)) + (\[Sigma]k*(2 + ph1)));
A14 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*btMw)/(\[Rho]k*btk)))*(1 -
ph2)) + (ph2*\[Rho]Mi*btMi)/(\[Rho]k*
btk))) + (ph3*\[Rho]Ag*btAg)/(\[Rho]k*btk);
A15 = ((1 -
ph3)*((((1 - ph1) + ((ph1*\[Rho]Mw*CPMw)/(\[Rho]k*CPk)))*(1 -
ph2)) + (ph2*\[Rho]Mi*CPMi)/(\[Rho]k*
CPk))) + (ph3*\[Rho]Ag*CPAg)/(\[Rho]k*CPk);
A16 = (\[Kappa]Ag + ((m - 1)*C16*B16*\[Kappa]k) - ((m - 1)*
ph3*(C16*B16*\[Kappa]k - \[Kappa]Ag)))/(\[Kappa]Ag + ((m - 1)*
C16*B16*\[Kappa]k) + (ph3*(C16*B16*\[Kappa]k - \[Kappa]Ag)));
B16 = (\[Kappa]Mi + ((m - 1)*C16*\[Kappa]k) - ((m - 1)*
ph2*((C16*\[Kappa]k) - \[Kappa]Mi)))/(\[Kappa]Mi + ((m - 1)*
C16*\[Kappa]k) + (ph2*((C16*\[Kappa]k) - \[Kappa]Mi)));
C16 = (\[Kappa]Mw + ((m - 1)*\[Kappa]k) - ((m - 1)*
ph1*(\[Kappa]k - \[Kappa]Mw)))/(\[Kappa]Mw + ((m -
1)*\[Kappa]k) + (ph1*(\[Kappa]k - \[Kappa]Mw)));
(*Equations*)

With[{soln = {1 -> 1}}, u101[y_] := Sin[la*x]*g1[y] /. soln;
u102[y_] := (m1[y] + Cos[2*la*x]*n1[y]) /. soln;
u111[y_] := Sin[la*x]*g2[y] /. soln;
u112[y_] := (m2[y] + Cos[2*la*x]*n2[y]) /. soln;
u201[y_] := Sin[la*x]*g3[y] /. soln;
u202[y_] := (m3[y] + Cos[2*la*x]*n3[y]) /. soln;
u211[y_] := Sin[la*x]*g4[y] /. soln;
u212[y_] := (m4[y] + Cos[2*la*x]*n4[y]) /. soln;
u10[y_] := (u100[y] + ep*u101[y] + ep^2*u102[y]) /. soln;
u11[y_] := (u110[y] + ep*u111[y] + ep^2*u112[y]) /. soln;
u20[y_] := (u200[y] + ep*u201[y] + ep^2*u202[y]) /. soln;
u21[y_] := (u210[y] + ep*u211[y] + ep^2*u212[y]) /. soln;
T101[y_] := Sin[la*x]*f1[y] /. soln;
T102[y_] := (h1[y] + Cos[2*la*x]*k1[y]) /. soln;
T111[y_] := Sin[la*x]*f2[y] /. soln;
T112[y_] := (h2[y] + Cos[2*la*x]*k2[y]) /. soln;
T201[y_] := Sin[la*x]*f3[y] /. soln;
T202[y_] := (h3[y] + Cos[2*la*x]*k3[y]) /. soln;
T211[y_] := Sin[la*x]*f4[y] /. soln;
T212[y_] := (h4[y] + Cos[2*la*x]*k4[y]) /. soln;
T10[y_] := (T100[y] + ep*T101[y] + ep^2*T102[y]) /. soln;
T11[y_] := (T110[y] + ep*T111[y] + ep^2*T112[y]) /. soln;
T20[y_] := (T200[y] + ep*T201[y] + ep^2*T202[y]) /. soln;
T21[y_] := (T210[y] + ep*T211[y] + ep^2*T212[y]) /. soln];
U[x_, y_] :=
Piecewise[{{(u20[y] + (\[Beta]*(E^(I*t))*u21[y])), \[Alpha] <= y <=
1}, {(u10[y] + (\[Beta]*(E^(I*t))*u11[y])), -1 <=
y <= \[Alpha]}}];

eq = {A1*(1 +
1/\[Beta]1)*(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*(1 + 1/\[Beta]1)*(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,
A11*(u200''[
y] + (u200'[y]/(y + k))) - ((A13*B13*C13)*(k/(y + k))^2*
M1^2)*u200[y] + ((A13*B13*C13*M1*b1) + (mu*l)) + (A14*Gr1*
T200[y]) == 0,
A11*(u210''[
y] + (u210'[y]/(y + k))) - (((A13*B13*C13)*(k/(y + k))^2*
M1^2) + (I*A12*H1^2))*
u210[y] + (mu*l1) + (A14*Gr1*T210[y]) == 0,
A1*(1 + 1/\[Beta]1)*(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*(1 + 1/\[Beta]1)*(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,
A11*(g3''[y] + (g3'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g3[y])) - ((A13*B13*C13)*(k/(y + k))^2*M1^2)*
g3[y] + (A14*Gr1*f3[y]) == 0,
A11*(g4''[y] + (g4'[y]/(y + k)) - (la^2*(k/(y + k))^2*
g4[y])) - (((A13*B13*C13)*(k/(y + k))^2*M1^2) + (I*A12*
H1^2))*g4[y] + (A14*Gr1*f4[y]) == 0,
A1*(1 + 1/\[Beta]1)*(m1''[
y] + (m1'[y]/(y + k))) - ((A3*B3*C3)*(k/(y + k))^2*M^2)*
m1[y] + (A4*Gr*h1[y]) == 0,
A1*(1 + 1/\[Beta]1)*(n1''[
y] + (n1'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*n1[y])) - ((A3*
B3*C3)*(k/(y + k))^2*M^2)*n1[y] + (A4*Gr*k1[y]) == 0,
A1*(1 + 1/\[Beta]1)*(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*(1 + 1/\[Beta]1)*(n2''[
y] + (n2'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n2[y])) - (((A3*B3*C3)*(k/(y + k))^2*M^2) + (I*A2*H^2))*
n2[y] + (A4*Gr*k2[y]) == 0,
A11*(m3''[y] + (m3'[y]/(y + k))) - ((A13*B13*C13)*(k/(y + k))^2*
M1^2)*m3[y] + (A14*Gr1*h3[y]) == 0,
A11*(n3''[y] + (n3'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n3[y])) - ((A13*B13*C13)*(k/(y + k))^2*M1^2)*
n3[y] + (A14*Gr1*k3[y]) == 0,
A11*(m4''[y] + (m4'[y]/(y + k))) - (((A13*B13*C13)*(k/(y + k))^2*
M1^2) + (I*A12*H1^2))*m4[y] + (A14*Gr1*h4[y]) == 0,
A11*(n4''[y] + (n4'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
n4[y])) - (((A13*B13*C13)*(k/(y + k))^2*M1^2) + (I*A12*
H1^2))*n4[y] + (A14*Gr1*k4[y]) ==
0, (((A6*B6*C6)/
Pr)*(T100''[y] + (T100'[y]/(y + k)))) + (\[CapitalOmega]*
T100[y]) ==
0, (((A6*B6*C6)/
Pr)*(T110''[
y] + (T110'[y]/(y + k)))) + ((\[CapitalOmega] - (I*A5*H^2))*
T110[y]) ==
0, (((A16*B16*C16)/
Pr1)*(T200''[y] + (T200'[y]/(y + k)))) + (\[CapitalOmega]1*
T200[y]) ==
0, (((A16*B16*C16)/
Pr1)*(T210''[
y] + (T200'[y]/(y + k)))) + ((\[CapitalOmega]1 - (I*A15*
H1^2))*T210[y]) ==
0, (((A6*B6*C6)/
Pr)*(f1''[y] + (f1'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f1[y]))) + (\[CapitalOmega]*f1[y]) ==
0, (((A6*B6*C6)/
Pr)*(f2''[y] + (f2'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f2[y]))) + ((\[CapitalOmega] - (I*A5*H^2))*f2[y]) ==
0, (((A16*B16*C16)/
Pr1)*(f3''[
y] + (f3'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f3[y]))) + (\[CapitalOmega]1*f3[y]) ==
0, (((A16*B16*C16)/
Pr1)*(f4''[
y] + (f4'[y]/(y + k)) - (la^2*(k/(y + k))^2*
f4[y]))) + ((\[CapitalOmega]1 - (I*A15*H1^2))*f4[y]) ==
0, (((A6*B6*C6)/
Pr)*(h1''[y] + (h1'[y]/(y + k)))) + (\[CapitalOmega]*h1[y]) ==
0, (((A6*B6*C6)/
Pr)*(k1''[y] + (k1'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k1[y]))) + (\[CapitalOmega]*k1[y]) ==
0, (((A6*B6*C6)/
Pr)*(h2''[y] + (h2'[y]/(y + k)))) + ((\[CapitalOmega] - (I*
A5*H^2))*h2[y]) ==
0, (((A6*B6*C6)/
Pr)*(k2''[y] + (k2'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k2[y]))) + ((\[CapitalOmega] - (I*A5*H^2))*k2[y]) ==
0, (((A16*B16*C16)/
Pr1)*(h3''[y] + (h3'[y]/(y + k)))) + (\[CapitalOmega]1*
h3[y]) ==
0, (((A16*B16*C16)/
Pr1)*(k3''[
y] + (k3'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k3[y]))) + (\[CapitalOmega]1*k3[y]) ==
0, (((A16*B16*C16)/
Pr1)*(h4''[
y] + (h4'[y]/(y + k)))) + ((\[CapitalOmega]1 - (I*A15*
H1^2))*h4[y]) ==
0, (((A16*B16*C16)/
Pr1)*(k4''[
y] + (k4'[y]/(y + k)) - (4*la^2*(k/(y + k))^2*
k4[y]))) + ((\[CapitalOmega]1 - (I*A15*H1^2))*k4[y]) ==
0};
(*Boundary conditions*)

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[\[Alpha]] == u200[\[Alpha]],
u110[\[Alpha]] ==
u210[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(u100'[\[Alpha]]) == (u200'[\[Alpha]]), (mu*(1 \
+ 1/\[Beta]1))*(u110'[\[Alpha]]) == (u210'[\[Alpha]]),
g1[\[Alpha]] == g3[\[Alpha]],
g2[\[Alpha]] ==
g4[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(g1'[\[Alpha]]) == (g3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(g2'[\[Alpha]]) == (g4'[\[Alpha]]),
m1[\[Alpha]] == m3[\[Alpha]], m2[\[Alpha]] == m4[\[Alpha]],
n1[\[Alpha]] == n3[\[Alpha]],
n2[\[Alpha]] ==
n4[\[Alpha]], (mu*(1 +
1/\[Beta]1))*(m1'[\[Alpha]]) == (m3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(m2'[\[Alpha]]) == (m4'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(n1'[\[Alpha]]) == (n3'[\[Alpha]]), (mu*(1 +
1/\[Beta]1))*(n2'[\[Alpha]]) == (n4'[\[Alpha]]),
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[\[Alpha]] == T200[\[Alpha]], f1[\[Alpha]] == f3[\[Alpha]],
h1[\[Alpha]] == h3[\[Alpha]], k1[\[Alpha]] == k3[\[Alpha]],
T110[\[Alpha]] == T210[\[Alpha]], f2[\[Alpha]] == f4[\[Alpha]],
h2[\[Alpha]] == h4[\[Alpha]],
k2[\[Alpha]] ==
k4[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
T100'[\[Alpha]] ==
T200'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
f1'[\[Alpha]] ==
f3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*h1'[\[Alpha]] ==
h3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
k1'[\[Alpha]] ==
k3'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
T110'[\[Alpha]] ==
T210'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
f2'[\[Alpha]] ==
f4'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*h2'[\[Alpha]] ==
h4'[\[Alpha]], ((\[Kappa]*A6*B6*C6)/(A16*B16*C16))*
k2'[\[Alpha]] == k4'[\[Alpha]]};

(*Solution Methodology*)

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"},
"StartingInitialConditions" ->
Table[var0[[i]] == 0 /. y -> -1, {i, Length[var]}]};
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,
Method -> {"Chasing", "ChasingType" -> "LinearChasing",
"ExtraPrecision" -> 10}, MaxSteps -> Infinity,
AccuracyGoal -> 10, PrecisionGoal -> 10];


Please pay attention that we use Method -> {"Chasing", "ChasingType" -> "LinearChasing","ExtraPrecision" -> 10} instead of shooting method. Visualization

{Plot[Evaluate[Re[var0] /. sol[[1]]], {y, -1, 1}, PlotRange -> All,
PlotLegends -> var, Frame -> True, FrameLabel -> Automatic],
p1 = Plot[U[0, y] /. sol[[1]] // Re, {y, -1, 1}, PlotRange -> All],
p2 = Plot[
Piecewise[{{Re[(T20[y] + (\[Beta]*(E^(I*t))*T21[y]))], \[Alpha] <=
y <= 1}, {Re[(T10[y] + (\[Beta]*(E^(I*t))*T11[y]))], -1 <=
y <= \[Alpha]}}] /. sol[[1]], {y, -1, 1}, PlotRange -> All]}


Update 1. Using Chasing method we can compute up to M^2=376 without warning from NDSolve and plot pictures as follows

While with shooting method we have in this case

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.

NDSolve::berr: The scaled boundary value residual error of 562.4155311689224 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.
`

Plotting results look like this

• sir, Here in the solution you are using chasing instead of shooting technique. What is the difference sir? In what ways it will effect solution Commented Oct 2, 2023 at 4:26
• Please, see Update 1 to my answer. Commented Oct 2, 2023 at 5:42
• okay sir, Thank you sir Commented Oct 2, 2023 at 7:09