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

Question

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/… — Alex Trounev, Commented Oct 1, 2023 at 6:26
@AlexTrounev Okay, sir, I will edit my accordingly. Thank you so much for the clarity sir — Komal Goyal, 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. — Alex Trounev, 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. — Komal Goyal, Commented Oct 1, 2023 at 8:36

Alex Trounev · Accepted Answer · 2023-10-02 05:40:31Z

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 — Komal Goyal, Commented Oct 2, 2023 at 4:26

Stack Exchange Network

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

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
differential-equations
equation-solving
boundary-conditions
or ask your own question.

Linked

Hot Network Questions

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

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged differential-equationsequation-solvingboundary-conditions or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
differential-equations
equation-solving
boundary-conditions
or ask your own question.