How to minimize numerical error in the code while solving the system of ODEs

I just want to ask in the graph of the code mentioned in the link, while plotting u from {-0.0001,y,0.0001}, the jump is coming because of the numerical error u1-u2. Can anyone please tell me, how to minimise that error? So that the graph will become continuous.

Solving a system of ODE using shooting technique with NDSolve command

Code

$Beta] = 10^-3; \[Epsilon] = 1/10; \[Alpha] = 1; \[Zeta] = 0; H = 1; M = 1; k = 3/2; z = \[Pi]/2; Subscript[\[Phi], 1] = 10^-2; Subscript[\[Phi], 2] = 10^-2; \[Rho]w = 997 + 1/10; \[Mu]w = 1004*10^-6; \[Sigma]w = 11/2*10^-6; \[Rho]k = 783; \[Mu]k = 151 10^-5; \[Sigma]k = 5*10^-11; \[Rho]fe = 5200; \[Sigma]fe = 25000; \[Rho]Mi = 5060; \[Sigma]Mi = 209*10^2; t = \[Pi]/3; Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^(5/2); 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])^(5/2)*(1 - Subscript[\[Phi], 1])^(5/ 2)); Subscript[B, 2] = (((1 - Subscript[\[Phi], 1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 - Subscript[\[Phi], 2])) + (Subscript[\[Phi], 2]*\[Rho]Mi)/\[Rho]k; Subscript[B, 3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w* Subscript[A, 3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 - Subscript[\[Phi], 2])) + (\[Sigma]w* Subscript[A, 3]*(2 + Subscript[\[Phi], 2]))); \[Sigma] = \[Sigma]w/\[Sigma]k; \[Mu] = \[Mu]w/\[Mu]k; \[Rho] = \[Rho]w/\[Rho]k; M1 = Sqrt[\[Mu]/\[Sigma]]*M; H1 = Sqrt[\[Mu]/\[Sigma]]*H; eq = {(Subscript[A, 1]/ Subscript[A, 2]*(u100''[ y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2* u100[y]))) - ((Subscript[A, 3]*M^2)/Subscript[A, 2]* u100[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(u110''[ y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2* u110[y]))) - (((Subscript[A, 3]*M^2)/Subscript[A, 2] + I*H^2)*u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(u200''[ y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2* u200[y]))) - ((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]* M1^2)/Subscript[B, 2]* u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(u210''[ y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2* u210[y]))) - (((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]* M1^2)/Subscript[B, 2] + I*H1^2)* u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/ Subscript[B, 2] + (I*H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/ Subscript[B, 2] + (I*H1^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g1''[y] + (1/(y + k)*g1'[y]))) - (g1[ y]*((Subscript[A, 1]/ Subscript[A, 2]*1/(y + k)^2) + (Subscript[A, 3]*M^2)/ Subscript[A, 2])) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[A, 3]* M^2)/Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[ y]*((Subscript[A, 1]/ Subscript[A, 2]*1/(y + k)^2) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[ y]*((Subscript[B, 1]/ Subscript[B, 2]*1/(y + k)^2) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0}; U[-1] = ((u200[ y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[ y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))) /. y -> -1; U[1] = ((u100[ y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[ y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. y -> 1; U1 = ((u200[ y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[ y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))) /. y -> 0; U2 = ((u100[ y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[ y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. y -> 0; bc = {U[1] == 0, u110[1] == 0, U[-1] == 0, u210[-1] == 0, f2[1] == 0, f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1], f5[-1] == -(Cos[\[Zeta]]*u200'[-1]), f6[-1] == -(Sin[\[Zeta]]*u200'[-1]), f7[-1] == -(Cos[\[Zeta]]*u210'[-1]), f8[-1] == -(Sin[\[Zeta]]*u210'[-1]), g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]), g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]), g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]), g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]), g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4* u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4* u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]), g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4* u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4* u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4* u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]), g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4* u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), U1 == U2, 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.2}; op = Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}}; var0 = Table[var[[i]]@y, {i, Length[var]}]; var1 = D[var0, y]; var2 = D[var1, y]; sol0 = Solve[eq, var2]; sol01 = sol0 /. y -> 1; sol02 = sol0 /. y -> -1; bc1 = bc /. sol01[[1]] /. sol02[[1]];  Numerical solution sol = NDSolve[{eq, bc1}, var, coord, op];  Boundary conditions test bc1 /. sol[[1]] Out[]= {False, False, False, True, False, False, False, False, \ True, True, True, False, False, False, False, False, False, False, \ True, True, True, True, True, True, False, False, False, False, \ False, False, False, False, False, False, False, False, False, False, \ False, False, False, False, False, False, False, False, False, False}  Velocity jump at y=0 (U1 - U2) /. sol[[1]] Out[]= 1.62451*10^-7 + 1.48518*10^-10 I  The question is can we reduce the real part jump on one order to 10^-8? Velocity jump visualization With[{n1 = sol[[1]]}, u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1; u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1; u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1; u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1; u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z]) /. n1; u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z]) /. n1; u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z]) /. n1; u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z]) /. n1; u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. n1; u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. n1; u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. n1; u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. n1]; p1 = Plot[ Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}]  Visualization in 2D • Please make your question self contained by including code and a picture Commented Jun 9, 2023 at 13:19 • Thank you, @MarcoB, for the suggestion, Actually the code is a little big, and already executed on Mathematica Stack Exchange. So, I think this is a better way to present the code. Kindly see by one simple click, please feel free to ask if you have any doubts, thank you so much. Commented Jun 9, 2023 at 16:59 • @KomalGoyal In this model jump is about 1.62451*10^-7 only. Why do you need 10^-8? Commented Jun 11, 2023 at 2:31 • @AlexTrounev, Because in the present case, if we plot 3-D function u(y,z) because of that jump that 3-D graph is discontinuous, but I want a continuous 3-D graph which will come when a numerical error is less than 10^-8 Commented Jun 11, 2023 at 3:45 • @KomalGoyal The best what I got with NDSolve is u1-u2=3.83437*10^-8 + 3.25555*10^-11 I. Commented Jun 11, 2023 at 8:40 1 Answer We can decrease error in 4.2 times by using NDSolve twice as follows \[Beta] = 10^-3; \[Epsilon] = 1/10; \[Alpha] = 1; \[Zeta] = 0; H = 1; M = 1; k = 3/2; Subscript[\[Phi], 1] = 10^-2; Subscript[\[Phi], 2] = 10^-2; \[Rho]w = 997 + 1/10; \[Mu]w = 1004*10^-6; \[Sigma]w = 11/2*10^-6; \[Rho]k = 783; \[Mu]k = 151 10^-5; \[Sigma]k = 5*10^-11; \[Rho]fe = 5200; \[Sigma]fe = 25000; \[Rho]Mi = 5060; \[Sigma]Mi = 209*10^2; t = \[Pi]/3; Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^(5/2); 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])^(5/2)*(1 - Subscript[\[Phi], 1])^(5/ 2)); Subscript[B, 2] = (((1 - Subscript[\[Phi], 1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 - Subscript[\[Phi], 2])) + (Subscript[\[Phi], 2]*\[Rho]Mi)/\[Rho]k; Subscript[B, 3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w* Subscript[A, 3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 - Subscript[\[Phi], 2])) + (\[Sigma]w* Subscript[A, 3]*(2 + Subscript[\[Phi], 2]))); \[Sigma] = \[Sigma]w/\[Sigma]k; \[Mu] = \[Mu]w/\[Mu]k; \[Rho] = \[Rho]w/\[Rho]k; M1 = Sqrt[\[Mu]/\[Sigma]]*M; H1 = Sqrt[\[Mu]/\[Sigma]]*H; eq = {(Subscript[A, 1]/ Subscript[A, 2]*(u100''[ y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2* u100[y]))) - ((Subscript[A, 3]*M^2)/Subscript[A, 2]* u100[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(u110''[ y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2* u110[y]))) - (((Subscript[A, 3]*M^2)/Subscript[A, 2] + I*H^2)*u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(u200''[ y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2* u200[y]))) - ((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]* M1^2)/Subscript[B, 2]* u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(u210''[ y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2* u210[y]))) - (((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]* M1^2)/Subscript[B, 2] + I*H1^2)* u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/ Subscript[B, 2] + (I*H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/ Subscript[B, 2] + (I*H1^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g1''[y] + (1/(y + k)*g1'[y]))) - (g1[ y]*((Subscript[A, 1]/ Subscript[A, 2]*1/(y + k)^2) + (Subscript[A, 3]*M^2)/ Subscript[A, 2])) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[A, 3]* M^2)/Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2]))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[ y]*((Subscript[A, 1]/ Subscript[A, 2]*1/(y + k)^2) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[A, 1]/ Subscript[A, 2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[ y]*((Subscript[A, 1]/ Subscript[A, 2]*(1/(y + k)^2 + 4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/ Subscript[A, 2] + (I*H^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2)) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[ y]*((Subscript[B, 1]/ Subscript[B, 2]*1/(y + k)^2) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0, (Subscript[B, 1]/ Subscript[B, 2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[ y]*((Subscript[B, 1]/ Subscript[B, 2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]* Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I* H1^2)))) == 0}; U[-1] = ((u200[ y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[ y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))) /. y -> -1; U[1] = ((u100[ y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[ y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. y -> 1; U1 = ((u200[ y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[ y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))) /. y -> 0; U2 = ((u100[ y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[ y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. y -> 0; bc = {U[1] == 0, u110[1] == 0, U[-1] == 0, u210[-1] == 0, f2[1] == 0, f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1], f5[-1] == -(Cos[\[Zeta]]*u200'[-1]), f6[-1] == -(Sin[\[Zeta]]*u200'[-1]), f7[-1] == -(Cos[\[Zeta]]*u210'[-1]), f8[-1] == -(Sin[\[Zeta]]*u210'[-1]), g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]), g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]), g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]), g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]), g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4* u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4* u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]), g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4* u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4* u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4* u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]), g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4* u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), U1 == U2, 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.2}; 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]];  First numerical solution sol = NDSolve[{eq, bc1}, var, coord, op];  With this solution we have (U1 - U2) /. sol[[1]] (*Out[]= 1.62451*10^-7 + 1.48518*10^-10 I*)  Second numerical solution ic0 = var0 /. sol[[1]] /. y -> -1; ic01 = var1 /. sol[[1]] /. y -> -1; ic1 = Join[Table[var0[[i]] == ic0[[i]] /. y -> -1, {i, Length[var0]}], Table[var1[[i]] == ic01[[i]] /. y -> -1, {i, Length[var1]}]]; sol1 = NDSolve[{eq, ic1}, var, coord, Method -> {"DiscontinuityProcessing" -> False, "TimeIntegration" -> "Extrapolation"}, StartingStepSize -> 2 10^-13];  With this solution we compute (U1 - U2) /. sol1[[1]] Out[]= 3.83437*10^-8 + 3.25555*10^-11 I  Visualization With[{n1 = sol[[1]]}, u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1; u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1; u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1; u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1; u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z]) /. n1; u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z]) /. n1; u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z]) /. n1; u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z]) /. n1; u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. n1; u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. n1; u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. n1; u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. n1]; p1 = Plot[ Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}, Frame -> True] With[{n1 = sol1[[1]]}, u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1; u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1; u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1; u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1; u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + g3[y]*Cos[2*\[Alpha]*z]) /. n1; u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z]) /. n1; u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + g9[y]*Cos[2*\[Alpha]*z]) /. n1; u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z]) /. n1; u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. n1; u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. n1; u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. n1; u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. n1]; p2 = Plot[ Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}, PlotStyle -> {Red, Dashed}] Show\[p1, p2$;][1]][1]


To plot velocity as a function y, z we use

f[y_, z_] :=
If[y < 0, ((u200[
y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] +
f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] +
g8[y]*Sin[2*\[Alpha]*z] +
g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] +
f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] +
g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))),
((u100[
y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] +
f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] +
g2[y]*Sin[2*\[Alpha]*z] +
g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] +
f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] +
g5[y]*Sin[2*\[Alpha]*z] +
g6[y]*Cos[2*\[Alpha]*z]))))] /.
NDSolve[{eq, bc1} /. z -> z, var, coord, op][[1]];

p2 = Plot[
Evaluate[Table[Re[f[y, z]], {z, 0, 2 Pi, Pi/16}]], {y, -1, 1},
PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"y", "U"},
ImageSize -> 260,
LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold},
PlotLegends -> Table[z, {z, 0, 2 Pi, Pi/16}]]


To extract data from p2 we use

Z = Table[z, {z, 0, 2 Pi, Pi/16}] // N;
lst = Table[{Z[[i]],
Extract[p2[[1, 1]][[1]][[1]][[2 + i]][[1]][[2]], 1]}, {i,
Length[Z]}];
YY = Table[lst[[i, 2]][[All, 1]], {i, Length[Z]}];
uu = Table[lst[[i, 2]][[All, 2]], {i, Length[Z]}];
dimY = Table[Length[YY[[i]]], {i, Length[Z]}];

u = Interpolation[
Flatten[Table[{{YY[[i, j]], Z[[i]]}, uu[[i, j]]}, {i,
Length[Z]}, {j, Min[dimY]}], 1], InterpolationOrder -> 1];

Plot3D[u[y, z], {y, -1, 1}, {z, 0, 2 Pi}, Mesh -> None,
ColorFunction -> "Aquamarine", PlotTheme -> "Marketing",
AxesLabel -> {"y", "z", "U"}]


Second method

Y1 = YY[[1]]; uz = ConstantArray[0, {Length[Y1], Length[Z]}];

Do[uz[[All, i]] =
Table[If[
y < 0, ((u200[
y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] +
f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] +
g8[y]*Sin[2*\[Alpha]*z] +
g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] +
f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] +
g11[y]*Sin[2*\[Alpha]*z] +
g12[y]*Cos[2*\[Alpha]*z])))),
((u100[
y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] +
f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] +
g2[y]*Sin[2*\[Alpha]*z] +
g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] +
f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] +
g5[y]*Sin[2*\[Alpha]*z] +
g6[y]*Cos[2*\[Alpha]*z]))))] /. z -> Z[[i]], {y,
Y1}] /. NDSolve[{eq, bc1} /. z -> Z[[i]], var, coord,
op][[1]];, {i, Length[Z]}]

u2 = Interpolation[
Flatten[Table[{{Y1[[j]], Z[[i]]}, uz[[j, i]]}, {i, Length[Z]}, {j,
dimY1[[2]]}], 1]]

Plot3D[u2[y, z] // Re, {y, -1, 1}, {z, 0, Pi}, Mesh -> None,
ColorFunction -> "Aquamarine", PlotTheme -> "Marketing",
AxesLabel -> {"y", "z", "U"}]


• Thank You sir for time and consideration @AlexTrounev Commented Jun 11, 2023 at 12:41
• @AlexTrounvev, can you please see the 3-D profiles for the same code, as mentioned above. Commented Jun 12, 2023 at 5:40
• @KomalGoyal Do you mean Plot3D using z as coordinate? Commented Jun 12, 2023 at 6:09
• @AlexTrounvev yes sir, I tried two times NDSolve also Commented Jun 12, 2023 at 7:11
• @KomalGoyal Could you stop editing my answer? Commented Jun 12, 2023 at 17:28