BVP4c for three boundary set of boundary conditions

Respected sir/maam,

I was trying to solve my coupled non-linear ODEs with three intervals of boundary conditions. Here It is error "uncognized function or variable D1" and some minor errors, which I am not able to rectify. Can you please help. If any information is required, please feel free to ask.

Thanking You in Advance

You help will be highly appreciable

xb = 0.4;
xc = 0.8;
xmesh = [linspace(0,0.4,100),linspace(0.4,0.8,100),linspace(0.8,1,100)];
yinit = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1 ];
solinit = bvpinit(xmesh,yinit);

%constant
r=1;
m=1;
cp=1;
k=1;
s=1;
bt=1;
bc2=1;
r1=1;
m1=1;
cp1=1;
k1=1;
s1=1;
bt1=1;
bc1=1;
t=1;
L=1;
e=0.001;
H=1;
R=0.2;
Da1=0.5;
Da2=0.5;
Da3=0.5;
B=1;
M=1;
Gr=1;
Gc=1;
th=pi/2;
d=1;
Pr=1;
Ec=1;
a1=1;
a2=1;
a3=1;
f1=1;
Nb=1;
Nt=1;
Le=1;
g=1;
K1=0.001;
Bi=1;
n=1;
Db=1;
Dt=1;
K=1;
Pr1=21;
v=(m*r1)/(m1*r);
si=(s1/s);
mu=(m1/m);
ro=(r/r1);
btk=(bt1/bt);
btc=(bc1/bc2);
CP=(cp/cp1);
roc=(r*cp)/(r1*cp1);
k2=(k/k1);
a=(r*cp*k1)/(k*r*cp1);

H1=sqrt(v)*H;

R1=v*R;

M1=sqrt(si/mu)*M;
Gr1=btk*v*Gr;
Gc1=btc*v*Gc;

Ec1=CP*Ec;

Nt1=(Nt)/(Dt*a);
Nb1=(Nb)/(Db*a);
Le1=a*Le*Db;
g1=(v*g)/K;
K2=(v*K1)/K;
D1=-(1/Da1)-(M^2)

D1 = -3

D2=sin(th)*Gr;
D3=sin(th)*Gc;
D4=d*((R/H)^2);
D5=(1/Pr);
D6=(f1-a1)/Pr;
D7=((Ec/Da1)+((M^2)*Ec));
D8=(Nb/Pr);
D9=(Nt/Pr);
D10=Le*Pr*(H^2);
D11=Le*Pr*R;
D12=(Nt/Nb);
D13=g*Le*Pr;
D14=K1*Le*Pr;
D15=ro*((H1)^2);
D16=-(1/Da2)-((M1)^2);
D17=sin(th)*Gr1;
D18=sin(th)*Gc1;
D19=(1/Pr1);
D20=d*(((R1)/(H1))^2);
D21=((f1-2)*k2)/Pr;
D22=((Ec1/Da2)+(Ec1*((M1)^2)));
D23=(ro*CP*Nb1)/(Pr1);
D24=(ro*CP*Nt1)/(Pr1);
D25=Le1*Pr1*(H1)^2;
D26=Le1*Pr1*R1;
D27=Nt1/Nb1;
D28=g1*Le1*Pr1;
D29=K2*Le1*Pr1;
D30=-(1/Da3)-((M)^2);
D31=(f1-a3)/Pr;
D32=(Ec/Da3)+((M^2)*Ec);

hold on

for R = 0.2
    p = [r m cp k s bt bc2 r1 m1 cp1 k1 s1 bt1 bc1 t L e H R Da1 Da2 Da3 B M Gr Gc th d Pr Ec a1 a2 a3 f1 Nb Nt Le g K1 Bi n Db Dt K Pr1 D1 D2 D3
  D4 D5 D6 D7 D8 D9 D10 D11 D12 D13 D14 D15 D16 D17 D18 D19 D20 D21 D22 D23 D24 D25 D26 D27 D28 D29 D30 D31 D32];
    sol = bvp5c(@(x,y,r) f(x,y,r,p), @(ya,yb)bc(ya,yb,p), solinit);
    plot(sol.x,real(sol.y(1,:)));
end

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

hold off

function dydx = f(x,y,region,p) % equations being solved
r = p(1);
m = p(2);
cp = p(3);
k = p(4);
s = p(5);
bt = p(6);
bc2 = p(7);
r1 = p(8);
m1 = p(9);
cp1 = p(10);
k1 = p(11);
s1 = p(12);
bt1 = p(13);
bc1 = p(14);
t = p(15);
L = p(16);
e = p(17);
H = p(18);
R = p(19);
Da1 = p(20);
Da2 = p(21);
Da3 = p(22);
B = p(23);
M = p(24);
Gr = p(25);
Gc = p(26);
th = p(27);
d = p(28);
Pr = p(29);
Ec = p(30);
a1 = p(31);
a2 = p(32);
a3 = p(33);
f1 = p(34);
Nb = p(35);
Nt = p(36);
Le = p(37);
g = p(38);
K1 = p(39);
Bi = p(40);
n = p(41);
Db = p(42);
Dt = p(43);
K = p(44);
Pr1 = p(45);
D1 = p(46);
  D2 = p(47);
  D3 = p(48);
  D4 = p(49);
  D5 = p(50);
  D6 = p(51);
  D7 = p(52);
  D8 = p(53);
  D9 = p(54);
  D10 = p(55);
  D11 = p(56);
  D12 = p(57);
  D13 = p(58);
  D14 = p(59);
  D15 = p(60);
  D16 = p(61);
  D17 = p(62);
  D18 = p(63);
  D19 = p(64);
  D20 = p(65);
  D21 = p(66);
  D22 = p(67);
  D23 = p(68);
  D24 = p(69);
  D25 = p(70);
  D26 = p(71);
  D27 = p(72);
  D28 = p(73);
  D29 = p(74);
  D30 = p(75);
  D31 = p(76);
  D32 = p(77);

dydx = zeros(6,1);
switch region 
    case 1     % x in [0,0.4]
        dydx = [y(2);
            (R*y(2))-(D1*y(1))-(D2*y(5))-(D3*y(9))-(H^2);
            y(4);
            (R*y(4))-((D1-(1i*(H^2)))*y(3))-(D2*y(7))-(D3*y(11))-(H^2);
            y(6);
            (1/(D5-D4))*((R*y(6))+(D6*y(5))-(Ec*(y(2)^2))-(D7*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)));
            y(8);
            (1/(D5-D4))*((R*y(8))+((D6+(1i*(H^2)))*y(7))-(2*Ec*y(2)*y(4))-(2*D7*y(1)*y(3))-(D8*((y(6)*y(12))+((y(10)*y(8)))))-(2*D9*y(6)*y*(8)))
            y(10);
            (D11*y(10))+(D13*y(9))-(D12*((1/(D5-D4))*((R*y(6))+(D6*y(5))-(Ec*(y(2)^2))-(D7*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)))))+(D14)
            y(12);
            (D11*y(12))+((D13+(1i*D10))*y(11))-(D12*((1/(D5-D4))*((R*y(8))+((D6+(1i*(H^2)))*y(7))-(2*Ec*y(2)*y(4))-(2*D7*y(1)*y(3))-(D8*((y(6)*y(12))+((y(10)*y(8)))))-(2*D9*y(6)*y*(8)))))];
    case 2     % x in [0.4,0.8]
        dydx = [y(2);
            (B/(B+1))*((R1*y(2))-(D16*y(1))-(D17*y(5))-(D18*y(9))-D15);
            y(4);
            (B/(B+1))*((R1*y(4))-((D16-(1i*(H1^2)))*y(3))-(D17*y(7))-(D18*y(11))-D15);
            y(6);
           (1/(D19-D20))*((R1*y(6))+(D21*y(5))-(Ec1*((B+1)/B)*(y(2)^2))-(D22*(y(1)^2))-(D23*y(10)*y(6))-(D24*(y(6)^2)));
           y(8);
           (1/(D19-D20))*((R1*y(8))+((D21+(1i*(H1^2)))*y(7))-(2*Ec1*((B+1)/B)*y(2)*y(4))-(2*D22*y(1)*y(3))-(D23*((y(10)*y(8))+((y(12)*y(6)))))-(2*D24*y(6)*y(8)));
           y(10);     
           (D26*y(10))+(D28*y(9))-(D27*((1/(D19-D20))*((R1*y(6))+(D21*y(5))-(Ec1*((B+1)/B)*(y(2)^2))-(D22*(y(1)^2))-(D23*y(10)*y(6))-(D24*(y(6)^2)))))+D29;
           y(12);
           (D26*y(12))+(D28*y(11))-(D27*((1/(D19-D20))*((R1*y(8))+((D21+(1i*(H1^2)))*y(7))-(2*Ec1*((B+1)/B)*y(2)*y(4))-(2*D22*y(1)*y(3))-(D23*((y(10)*y(8))+((y(12)*y(6)))))-(2*D24*y(6)*y(8)))))
           ];
      case 3    % x in [0.8,1]
        dydx = [y(2);
            (R*y(2))-(D30*y(1))-(D2*y(5))-(D3*y(9))-(H^2);
            y(4);
            (R*y(4))-((D30-(1i*(H^2)))*y(3))-(D2*y(7))-(D3*y(11))-(H^2);
            y(6);
            (1/(D5-D4))*((R*y(6))+(D31*y(5))-(Ec*(y(2)^2))-(D32*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)));
            y(8);
            (1/(D5-D4))*((R*y(8))+((D31+(1i*(H^2)))*y(7))-(2*Ec*y(4)*y(2))-(2*D32*y(1)*y(3))-(D8*((y(10)*y(8))+((y(6)*y(12)))))-(2*D9*y(6)*y*(8)));
            y(10);
            (D11*y(10))+(D13*y(9))-(D12*((1/(D5-D4))*((R*y(6))+(D31*y(5))-(Ec*(y(2)^2))-(D32*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)))))+(D14);
            y(12);
            (D11*y(12))+((D13+(1i*D10))*y(11))-(D12*((1/(D5-D4))*((R*y(8))+((D31+(1i*(H^2)))*y(7))-(2*Ec*y(4)*y(2))-(2*D32*y(1)*y(3))-(D8*((y(10)*y(8))+((y(6)*y(12)))))-(2*D9*y(6)*y*(8)))))
            ];
end     
end

function res = bc(YL,YR,p)
B = p(23);
m = p(9)/p(2);
k2 = p(4)/p(11);
Db = p(42);
res = [YL(1,1)            % u10(0)=0
    YL(3,1)            % u11(0)=0
    YL(5,1)            % T10(0)=0
    YL(7,1)            % T11(0)=0
    YL(9,1)            % Phi(0)=0
    YL(11,1)           % Phi(0)=0
    YR(1,1)-YL(1,2)    % u20(0.4)=u10(0.4)
    YR(3,1)-YL(3,2)    % u21(0.4)=u11(0.4)
    YR(1,2)-YL(1,3)    % u30(0.8)=u20(0.8)
    YR(3,2)-YL(3,3)    % u31(0.8)=u21(0.8)
    YR(5,1)-YL(5,2)    % T20(0.4)=T10(0.4)
    YR(7,1)-YL(7,2)    % T21(0.4)=T11(0.4)
    YR(5,2)-YL(5,3)    % T30(0.8)=T20(0.8)
    YR(7,2)-YL(7,3)    % T31(0.8)=T21(0.8)
    YR(9,1)-YL(9,2)    % P20(0.4)=P10(0.4)
    YR(11,1)-YL(11,2)  % P21(0.4)=P11(0.4)
    YR(9,2)-YL(9,3)    % P30(0.8)=P20(0.8)
    YR(11,2)-YL(11,3)              % P31(0.8)=P21(0.8)
    (m*((B+1)/B))*YR(2,1)-YL(2,2)  % m*(1+(1/b))*u20'(0.4)=u10'(0.4)
    (m*((B+1)/B))*YR(4,1)-YL(4,2)  % m*(1+(1/b))*u21'(0.4)=u11'(0.4)
    YR(2,2)-((m*((B+1)/B))*YL(2,3)) % m*(1+(1/b))*u20'(0.8)=u30'(0.8)
    YR(4,2)-((m*((B+1)/B))*YL(4,3))  % m*(1+(1/b))*u21'(0.8)=u31'(0.8)
    YR(6,1)-(k2*YL(6,2))             % T20'[0.4]=k*T10'[0.4]
    YR(8,1)-(k2*YL(8,2))             % T21'[0.4]=k*T11'[0.4]
    (k2*YR(6,2))-(YL(6,3))           % k*T30'[0.8]=T20'[0.8]
    (k2*YR(8,2))-(YL(8,3))           % k*T31'[0.8]=T21'[0.8]
    YR(10,1)-(Db*YL(10,2))           % P20'[0.4]=Db*P10'[0.4]
    YR(12,1)-(Db*YL(12,2))           % P21'[0.4]=Db*P11'[0.4]
    (Db*YR(10,2))-(YL(10,3))         % Db*P30'[0.8]=P20'[0.8]
    (Db*YR(12,2))-(YL(12,3))           % Db*T31'[0.8]=P21'[0.8]
    YR(1,3)            % u30(1)=0
    YR(3,3)            % u31(1)=0
    YR(5,3)            % T30(1)=0
    YR(7,3)            % T31(1)=0
    YR(9,3)            % Phi30(1)=0
    YR(11,3)           % Phi31(1)=0
    ];
end 

Answer 1

[Torsten]

You didn't try to understand the code I provided for the similar problem, did you ?

You still supply 6 initial conditions although you have to supply 12.

You still add the equation numbers for the different regions although they have to start with 1 in every region and go up until 12. Your method gives you a numbering up to 12 + 12 + 12 = 36 (1-12, 13-24, 25-36) for the equation numbers within the regions, although the index should only run from 1-12 in every region (1-12, 1-12, 1-12).

[Komal Goyal]

@Torsten sir, I was trying sir to undrstand that code, and write the similiar code. Thank you for your suggestion, I will edit as per your comments and try to get the output

Thanking You

[Komal Goyal]

@Torsten sir, I am trying my level best to learn this bvp4c coding, I did the changes in the ODEs, but unable to understand what are the changes required for the boundary conditions and where. Can you please suggest that once.

Thanking You

[Torsten]

The initial value vector is still a vector of length 6 instead of 12:

yinit = [1; 1; 1; 1; 1; 1 ];

All the D1,D2,...,D32 in function "f" are undefined.

[Komal Goyal]

@Torsten sir, I am trying to define D1-D32 in function f, but it is showing error regarding vertcat. Can you please suggest how to rectify this error.

[Torsten]

It works now, but it takes too long to run it here online.

R1, H1 and Ec1 were not defined. I set them as

R1 = R;
H1 = H;
Ec1 = Ec;

xb = 0.4;
xc = 0.8;
xmesh = [linspace(0,0.4,100),linspace(0.4,0.8,100),linspace(0.8,1,100)];
yinit = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1 ];
solinit = bvpinit(xmesh,yinit);

%constant
r=1;
m=1;
cp=1;
k=1;
s=1;
bt=1;
bc2=1;
r1=1;
m1=1;
cp1=1;
k1=1;
s1=1;
bt1=1;
bc1=1;
t=1;
L=1;
e=0.001;
H=1;
R=0.2;
Da1=0.5;
Da2=0.5;
Da3=0.5;
B=1;
M=1;
Gr=1;
Gc=1;
th=pi/2;
d=1;
Pr=1;
Ec=1;
a1=1;
a2=1;
a3=1;
f1=1;
Nb=1;
Nt=1;
Le=1;
g=1;
K1=0.001;
Bi=1;
n=1;
Db=1;
Dt=1;
K=1;
Pr1=21;
v=(m*r1)/(m1*r);
si=(s1/s);
mu=(m1/m);
ro=(r/r1);
btk=(bt1/bt);
btc=(bc1/bc2);
CP=(cp/cp1);
roc=(r*cp)/(r1*cp1);
k2=(k/k1);
a=(r*cp*k1)/(k*r*cp1);
H1=sqrt(v)*H;
R1=v*H;
M1=sqrt(si/mu)*M;
Gr1=btk*v*Gr;
Gc1=btc*v*Gc;
Ec1=CP*Ec;
Nt1=(Nt)/(Dt*a);
Nb1=(Nb)/(Db*a);
Le1=a*Le*Db;
g1=(v*g)/K;
K2=(v*K1)/K;
D1=-(1/Da1)-(M^2)
D2=sin(th)*Gr;
D3=sin(th)*Gc;
D4=d*((R/H)^2);
D5=(1/Pr);
D6=(f1-a1)/Pr;
D7=((Ec/Da1)+((M^2)*Ec));
D8=(Nb/Pr);
D9=(Nt/Pr);
D10=Le*Pr*(H^2);
D11=Le*Pr*R;
D12=(Nt/Nb);
D13=g*Le*Pr;
D14=K1*Le*Pr;
D15=ro*((H1)^2);
D16=-(1/Da2)-((M1)^2);
D17=sin(th)*Gr1;
D18=sin(th)*Gc1;
D19=(1/Pr1);
D20=d*(((R1)/(H1))^2);
D21=((f1-2)*k2)/Pr;
D22=((Ec1/Da2)+(Ec1*((M1)^2)));
D23=(ro*CP*Nb1)/(Pr1);
D24=(ro*CP*Nt1)/(Pr1);
D25=Le1*Pr1*(H1)^2;
D26=Le1*Pr1*R1;
D27=Nt1/Nb1;
D28=g1*Le1*Pr1;
D29=K2*Le1*Pr1;
D30=-(1/Da3)-((M)^2);
D31=(f1-a3)/Pr;
D32=(Ec/Da3)+((M^2)*Ec);

hold on
for R = 0.2
    p = [r m cp k s bt bc2 r1 m1 cp1 k1 s1 bt1 bc1 t L e H R Da1 Da2 Da3 B M Gr Gc th d Pr Ec a1 a2 a3 f1 Nb Nt Le g K1 Bi n Db Dt K Pr1 D1 D2 D3...
         D4 D5 D6 D7 D8 D9 D10 D11 D12 D13 D14 D15 D16 D17 D18 D19 D20 D21 D22 D23 D24 D25 D26 D27 D28 D29 D30 D31 D32];
    sol = bvp5c(@(x,y,r) f(x,y,r,p), @(ya,yb)bc(ya,yb,p), solinit);
    plot(sol.x,real(sol.y(1,:)));
end
hold off

function dydx = f(x,y,region,p) % equations being solved
r = p(1);
m = p(2);
cp = p(3);
k = p(4);
s = p(5);
bt = p(6);
bc2 = p(7);
r1 = p(8);
m1 = p(9);
cp1 = p(10);
k1 = p(11);
s1 = p(12);
bt1 = p(13);
bc1 = p(14);
t = p(15);
L = p(16);
e = p(17);
H = p(18);
R = p(19);
Da1 = p(20);
Da2 = p(21);
Da3 = p(22);
B = p(23);
M = p(24);
Gr = p(25);
Gc = p(26);
th = p(27);
d = p(28);
Pr = p(29);
Ec = p(30);
a1 = p(31);
a2 = p(32);
a3 = p(33);
f1 = p(34);
Nb = p(35);
Nt = p(36);
Le = p(37);
g = p(38);
K1 = p(39);
Bi = p(40);
n = p(41);
Db = p(42);
Dt = p(43);
K = p(44);
Pr1 = p(45);
D1 = p(46);
  D2 = p(47);
  D3 = p(48);
  D4 = p(49);
  D5 = p(50);
  D6 = p(51);
  D7 = p(52);
  D8 = p(53);
  D9 = p(54);
  D10 = p(55);
  D11 = p(56);
  D12 = p(57);
  D13 = p(58);
  D14 = p(59);
  D15 = p(60);
  D16 = p(61);
  D17 = p(62);
  D18 = p(63);
  D19 = p(64);
  D20 = p(65);
  D21 = p(66);
  D22 = p(67);
  D23 = p(68);
  D24 = p(69);
  D25 = p(70);
  D26 = p(71);
  D27 = p(72);
  D28 = p(73);
  D29 = p(74);
  D30 = p(75);
  D31 = p(76);
  D32 = p(77);

  R1 = R;
  H1 = H;
  Ec1 = Ec;
dydx = zeros(12,1);
switch region 
    case 1     % x in [0,0.4]
        dydx = [y(2);
               (R*y(2))-(D1*y(1))-(D2*y(5))-(D3*y(9))-(H^2);
               y(4);
               (R*y(4))-((D1-(1i*(H^2)))*y(3))-(D2*y(7))-(D3*y(11))-(H^2);
               y(6);
               (1/(D5-D4))*((R*y(6))+(D6*y(5))-(Ec*(y(2)^2))-(D7*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)));
               y(8);
               (1/(D5-D4))*((R*y(8))+((D6+(1i*(H^2)))*y(7))-(2*Ec*y(2)*y(4))-(2*D7*y(1)*y(3))-(D8*((y(6)*y(12))+((y(10)*y(8)))))-(2*D9*y(6)*y(8)));
               y(10);
               (D11*y(10))+(D13*y(9))-(D12*((1/(D5-D4))*((R*y(6))+(D6*y(5))-(Ec*(y(2)^2))-(D7*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)))))+(D14);
               y(12);
               (D11*y(12))+((D13+(1i*D10))*y(11))-(D12*((1/(D5-D4))*((R*y(8))+((D6+(1i*(H^2)))*y(7))-(2*Ec*y(2)*y(4))-(2*D7*y(1)*y(3))-(D8*((y(6)*y(12))+((y(10)*y(8)))))-(2*D9*y(6)*y(8)))))]
    case 2     % x in [0.4,0.8]
        dydx = [y(2);
               (B/(B+1))*((R1*y(2))-(D16*y(1))-(D17*y(5))-(D18*y(9))-D15);
               y(4);
               (B/(B+1))*((R1*y(4))-((D16-(1i*(H1^2)))*y(3))-(D17*y(7))-(D18*y(11))-D15);
               y(6);
               (1/(D19-D20))*((R1*y(6))+(D21*y(5))-(Ec1*((B+1)/B)*(y(2)^2))-(D22*(y(1)^2))-(D23*y(10)*y(6))-(D24*(y(6)^2)));
               y(8);
               (1/(D19-D20))*((R1*y(8))+((D21+(1i*(H1^2)))*y(7))-(2*Ec1*((B+1)/B)*y(2)*y(4))-(2*D22*y(1)*y(3))-(D23*((y(10)*y(8))+((y(12)*y(6)))))-(2*D24*y(6)*y(8)));
               y(10);     
               (D26*y(10))+(D28*y(9))-(D27*((1/(D19-D20))*((R1*y(6))+(D21*y(5))-(Ec1*((B+1)/B)*(y(2)^2))-(D22*(y(1)^2))-(D23*y(10)*y(6))-(D24*(y(6)^2)))))+D29;
               y(12);
               (D26*y(12))+(D28*y(11))-(D27*((1/(D19-D20))*((R1*y(8))+((D21+(1i*(H1^2)))*y(7))-(2*Ec1*((B+1)/B)*y(2)*y(4))-(2*D22*y(1)*y(3))-(D23*((y(10)*y(8))+((y(12)*y(6)))))-(2*D24*y(6)*y(8)))))];
      case 3    % x in [0.8,1]
        dydx = [y(2);
               (R*y(2))-(D30*y(1))-(D2*y(5))-(D3*y(9))-(H^2);
               y(4);
               (R*y(4))-((D30-(1i*(H^2)))*y(3))-(D2*y(7))-(D3*y(11))-(H^2);
               y(6);
               (1/(D5-D4))*((R*y(6))+(D31*y(5))-(Ec*(y(2)^2))-(D32*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)));
               y(8);
               (1/(D5-D4))*((R*y(8))+((D31+(1i*(H^2)))*y(7))-(2*Ec*y(4)*y(2))-(2*D32*y(1)*y(3))-(D8*((y(10)*y(8))+((y(6)*y(12)))))-(2*D9*y(6)*y(8)));
               y(10);
               (D11*y(10))+(D13*y(9))-(D12*((1/(D5-D4))*((R*y(6))+(D31*y(5))-(Ec*(y(2)^2))-(D32*(y(1)^2))-(D8*y(6)*y(10))-(D9*(y(6)^2)))))+(D14);
               y(12);
               (D11*y(12))+((D13+(1i*D10))*y(11))-(D12*((1/(D5-D4))*((R*y(8))+((D31+(1i*(H^2)))*y(7))-(2*Ec*y(4)*y(2))-(2*D32*y(1)*y(3))-(D8*((y(10)*y(8))+((y(6)*y(12)))))-(2*D9*y(6)*y(8)))))];
end
end

function res = bc(YL,YR,p)
B = p(23);
m = p(9)/p(2);
k2 = p(4)/p(11);
Db = p(42);
res = [YL(1,1)            % u10(0)=0
    YL(3,1)            % u11(0)=0
    YL(5,1)            % T10(0)=0
    YL(7,1)            % T11(0)=0
    YL(9,1)            % Phi(0)=0
    YL(11,1)           % Phi(0)=0
    YR(1,1)-YL(1,2)    % u20(0.4)=u10(0.4)
    YR(3,1)-YL(3,2)    % u21(0.4)=u11(0.4)
    YR(1,2)-YL(1,3)    % u30(0.8)=u20(0.8)
    YR(3,2)-YL(3,3)    % u31(0.8)=u21(0.8)
    YR(5,1)-YL(5,2)    % T20(0.4)=T10(0.4)
    YR(7,1)-YL(7,2)    % T21(0.4)=T11(0.4)
    YR(5,2)-YL(5,3)    % T30(0.8)=T20(0.8)
    YR(7,2)-YL(7,3)    % T31(0.8)=T21(0.8)
    YR(9,1)-YL(9,2)    % P20(0.4)=P10(0.4)
    YR(11,1)-YL(11,2)  % P21(0.4)=P11(0.4)
    YR(9,2)-YL(9,3)    % P30(0.8)=P20(0.8)
    YR(11,2)-YL(11,3)              % P31(0.8)=P21(0.8)
    (m*((B+1)/B))*YR(2,1)-YL(2,2)  % m*(1+(1/b))*u20'(0.4)=u10'(0.4)
    (m*((B+1)/B))*YR(4,1)-YL(4,2)  % m*(1+(1/b))*u21'(0.4)=u11'(0.4)
    YR(2,2)-((m*((B+1)/B))*YL(2,3)) % m*(1+(1/b))*u20'(0.8)=u30'(0.8)
    YR(4,2)-((m*((B+1)/B))*YL(4,3))  % m*(1+(1/b))*u21'(0.8)=u31'(0.8)
    YR(6,1)-(k2*YL(6,2))             % T20'[0.4]=k*T10'[0.4]
    YR(8,1)-(k2*YL(8,2))             % T21'[0.4]=k*T11'[0.4]
    (k2*YR(6,2))-(YL(6,3))           % k*T30'[0.8]=T20'[0.8]
    (k2*YR(8,2))-(YL(8,3))           % k*T31'[0.8]=T21'[0.8]
    YR(10,1)-(Db*YL(10,2))           % P20'[0.4]=Db*P10'[0.4]
    YR(12,1)-(Db*YL(12,2))           % P21'[0.4]=Db*P11'[0.4]
    (Db*YR(10,2))-(YL(10,3))         % Db*P30'[0.8]=P20'[0.8]
    (Db*YR(12,2))-(YL(12,3))           % Db*T31'[0.8]=P21'[0.8]
    YR(1,3)            % u30(1)=0
    YR(3,3)            % u31(1)=0
    YR(5,3)            % T30(1)=0
    YR(7,3)            % T31(1)=0
    YR(9,3)            % Phi30(1)=0
    YR(11,3)           % Phi31(1)=0
    ];
end 

[Komal Goyal]

Thank you so much sir for you help, fron the depth of my heart. Just last question I want to ask

If I want to plot y(1,:)+e^(1i*t)y(3,:) from first set of ODE in interval [0,0.4], y(1,:)+e^(1i*t)y(3,:) from second set of ODE in interval [0.4,0.8] and y(1,:)+e^(1i*t)y(3,:) from third set of ODE in interval [0.5,1], considering x-axis [0,1]. Hou can I do that?

[Torsten]

It's already done here:

https://de.mathworks.com/matlabcentral/answers/1904825-solving-system-of-coupled-boundary-value-problem-using-bvp4c?s_tid=srchtitle

plot(sol.x,real(sol.y(1,:)+exp(1i*t)*sol.y(3,:)));

[Komal Goyal]

@Torsten, Okay sir get it. Sir, R1, Ec1, H1 is already defined in list of constants. Why it is showing undefined sir? You can check my constant list, Just how I highlight R1, Ec1, H1 to make it clear.

[Torsten]

I don't know what you mean.

As I said, I added the three lines

R1 = R;
H1 = H;
Ec1 = Ec;

to your code to make it work.

I see the three constants now in the definitions in the calling program, but you didn't include them in the parameter list p for function f.

[Komal Goyal]

@Torsten, when I am exwcutiong the code, graph is coming fine, but this is the error which is coming at the last. How can we meet to the tolerance sir? because maximum error is very high

Warning: Unable to meet the tolerance without using more than 833 mesh points.

The last mesh of 300 points and the solution are available in the output argument.

The maximum error is 9.77749, while requested accuracy is 0.001.

> In bvp5c (line 279)

In final (line 111)

[Torsten]

You encounter the main problem in all numerical computations: the solver has problems getting a solution.

Either there are error in the problem formulation.

Or you made mistakes in coding your problem.

Or the problem is that hard the solver is unable to get a solution.

Or the problem does not have a solution.

To proceed further, you need knowledge about the solutions of your system that you expect.

Ask yourself:

Are the results I got so far plausible ? How can I improve the initial guesses for the solutions ?

I cannot help you in this respect because I don't know the underlying problem and thus cannot answer any of these questions.

[Komal Goyal]

Thank You @Torsten sir, I will work on that. I have one doubt that if we want to find the value if y(1,:) for a specific case out of the 3 cases so how to print that value sir? I am using fprintf(y(1,1)); but i think it is displaying the total of y(1,:) in all three regions. But I want for only first region.

Thanking you sir for the help in advance.

[Torsten]

y(1,1) is the value of y(1) in point xmesh(1), thus a scalar. Why do you write it is displaying the total of y(1,:) in all three regions ?

What exactly do you want to get ? The solution of y(1) in the first region ?

Then maybe

idx = sol.x<=0.4;
plot(sol.x(idx),real(sol.y(1,idx)))

?

[Komal Goyal]

Thank you so much sir for your help, yes I want exactly the same

idx = sol.x<=0.4;

plot(sol.x(idx),real(sol.y(1,idx)))

But if I want to find the value of the above expression at say x=1, then how to print that value sir?

[Torsten]

x = 1 does not exist in the above expression. You wanted to restrict y to the first region, and this region ends at x = 0.4.

Use "deval" to evaluate the solution at certain x-values:

x = 1;
y = deval(sol,x)

[Komal Goyal]

Respected @Torsten sir, I dont know, whether I am able to tell you propely or not what is my doubt. Here I attached my doubt that in which region I want to find values and code which I prepare because of your help. Kindly help me in that as well.