graph is the not plotted as desired because of tolerance level

hi,

I have a coupled non-linear differential equations

(d^2 f)/(dy^2 )+m2*g2*dB/dy-2*i*R2*g1*f - g3*G1*y - R4*g1 = 0

(d^2 B)/(dy^2 )+t4/(1-i*H1)*df/dy=0

(d^2 T)/(dy^2 )-1/2*g4*G1*PR*(f+ ̅f)+ER*PR*[g5*(df/dy*d ̅f/dy)+g6*m2*(dB/dy*dB̅/dy)]=0, where ̅f is conjugate of f and B̅ is conjugate of B

Boundary conditions are

f=0 at y=0

f=C1 at y=1

And

dB/dy-(t4/(P1* (1-i*H1 ) ))* B=0 at y=0

dB/dy+(t4/(P2 (1-i*H1 ) ))* B=0 at y=1

and

T=0 at y=0

T=1 at y=1

when i run the program, i get an error "Warning: Unable to meet the tolerance without using more than 1666 mesh

points. The last mesh of 833 points and the solution are available in the output argument. The maximum residual is 0.875363, while requested accuracy is 0.001. Unable to meet the tolerance without using more than 1666 mesh points", I tried using NMax but could not get the solution"

I have enclosed the program and graph for your reference, I need to get the graph similar to the one i have enclosed. Please help me in completing the graph.

Thank you

Mat lab program

close all
clc
p=0.1;
P1=2;
P2=2;
b1=0.00021;
b2=0.000058;
S1=0.005;
S2=580000;
G1=2;
EM2=20;
R1=997.1;
R2=5;
C1=1;
R3=4420;
H1=0.25;
K1=1;
R4=1;
PR=7.0; 
ER= 2.0;
cf=4179;
cs=0.56;
K2=0.613;
K3=7.2;
t1=(1./((1-p).^2.5));
t2=(1-p)+(p.*(R3./R1));
t3=(1-p)+p.*((R3.*b2)./(R1.*b1));
S=(S2./S1);
t4=1-((3*(1-S).*p)./((2+S)+(1-S).*p));
t5=(1-p)+(p.*R3.*cs)./(R1.*cf);
t6=(1+2.*(K2./K3)+2.*p.*(1-K2./K3))./(1+2.*(K2./K3)-p.*(1-K2./K3));
g1=t2./t1;
g2=1/t1;
g3=t3./t1;
g4=t5./t6;
g5=t1./t6;
g6=1./(t4.*t6);
m1=(t4./(P1.*(1-1i.*H1)));
m2=(t4./(P2.*(1-1i.*H1)));
dydx=@(x,y)[y(4);
    y(5);
    y(6);
    -EM2.*g2.*y(4)+2.*1i.*R2.*g1.*y(1)+g3.*G1.*x+R4.*g1;
    (-t4./(1-1i.*H1)).*y(3);
    1/2.*g4.*G1.*PR.*(y(1)+conj(y(1)))-ER.*PR.*(g5.*(y(4).*conj(y(4))+g6.*EM2.*(y(5).*conj(y(5)))))];
BC = @(ya,yb)[ya(1)-0;yb(1)-C1;ya(3)-0;yb(3)-1.0;ya(5)-m1.*ya(2);yb(5)+m2.*yb(2)];
yinit = [0.01;0.01;0.01;0.01;0.01;0.01];
solinit = bvpinit(linspace(0,1,50),yinit);
options = bvpset('AbsTol',1e-6);
U1 = bvp4c(dydx,BC,solinit,options);       
hold on
p=0.001;
P1=2;
P2=2;
b1=0.00021;
b2=0.000058;
S1=0.005;
S2=580000;
G1=2;
EM2=20;
R1=997.1;
R2=5;
C1=1;
R3=4420;
H1=0.5;
K1=1;
R4=1;
PR=7.0; 
ER= 2.0;
cf=4179;
cs=0.56;
K2=0.613;
K3=7.2;
t1=(1./((1-p).^2.5));
t2=(1-p)+(p.*(R3./R1));
t3=(1-p)+p.*((R3.*b2)./(R1.*b1));
S=(S2./S1);
t4=1-((3*(1-S).*p)./((2+S)+(1-S).*p));
t5=(1-p)+(p.*R3.*cs)./(R1.*cf);
t6=(1+2.*(K2./K3)+2.*p.*(1-K2./K3))./(1+2.*(K2./K3)-p.*(1-K2./K3));
g1=t2./t1;
g2=1/t1;
g3=t3./t1;
g4=t5./t6;
g5=t1./t6;
g6=1./(t4.*t6);
m1=(t4./(P1.*(1-1i.*H1)));
m2=(t4./(P2.*(1-1i.*H1)));
dydx=@(x,y)[y(4);
    y(5);
    y(6);
    -EM2.*g2.*y(4)+2.*1i.*R2.*g1.*y(1)+g3.*G1.*x+R4.*g1;
    (-t4./(1-1i.*H1)).*y(3);
    1/2.*g4.*G1.*PR.*(y(1)+conj(y(1)))-ER.*PR.*(g5.*(y(4).*conj(y(4))+g6.*EM2.*(y(5).*conj(y(5)))))];
BC = @(ya,yb)[ya(1)-0;yb(1)-C1;ya(3)-0;yb(3)-1.0;ya(5)-m1.*ya(2);yb(5)+m2.*yb(2)];
yinit = [0.01;0.01;0.01;0.01;0.01;0.01];
solinit = bvpinit(linspace(0,1,50),yinit);
options = bvpset('AbsTol',1e-6);
U2 = bvp4c(dydx,BC,solinit,options);       
hold on
p=0.02;
P1=2;
P2=2;
b1=0.00021;
b2=0.000058;
S1=0.005;
S2=580000;
G1=2;
EM2=20;
R1=997.1;
R2=5;
C1=1;
R3=4420;
H1=0.75;
K1=1;
R4=1;
PR=7.0; 
ER= 2.0;
cf=4179;
cs=0.56;
K2=0.613;
K3=7.2;
t1=(1./((1-p).^2.5));
t2=(1-p)+(p.*(R3./R1));
t3=(1-p)+p.*((R3.*b2)./(R1.*b1));
S=(S2./S1);
t4=1-((3*(1-S).*p)./((2+S)+(1-S).*p));
t5=(1-p)+(p.*R3.*cs)./(R1.*cf);
t6=(1+2.*(K2./K3)+2.*p.*(1-K2./K3))./(1+2.*(K2./K3)-p.*(1-K2./K3));
g1=t2./t1;
g2=1/t1;
g3=t3./t1;
g4=t5./t6;
g5=t1./t6;
g6=1./(t4.*t6);
m1=(t4./(P1.*(1-1i.*H1)));
m2=(t4./(P2.*(1-1i.*H1)));
dydx=@(x,y)[y(4);
    y(5);
    y(6);
    -EM2.*g2.*y(4)+2.*1i.*R2.*g1.*y(1)+g3.*G1.*x+R4.*g1;
    (-t4./(1-1i.*H1)).*y(3);
    1/2.*g4.*G1.*PR.*(y(1)+conj(y(1)))-ER.*PR.*(g5.*(y(4).*conj(y(4))+g6.*EM2.*(y(5).*conj(y(5)))))];
BC = @(ya,yb)[ya(1)-0;yb(1)-C1;ya(3)-0;yb(3)-1.0;ya(5)-m1.*ya(2);yb(5)+m2.*yb(2)];
yinit = [0.01;0.01;0.01;0.01;0.01;0.01];
solinit = bvpinit(linspace(0,1,50),yinit);
options = bvpset('AbsTol',1e-6);
U3 = bvp4c(dydx,BC,solinit,options);       
hold on
plot(U1.x,real([U1.y(3,:)]),'b','linewidth',1.5)
plot(U2.x,real([U2.y(3,:)]),'r','linewidth',1.5)
plot(U3.x,real([U3.y(3,:)]),'g','linewidth',1.5)

Answer 1

[Torsten]

dydx(5) = (-t4./(1-1i.*H1)).*y(4);

instead of

dydx(5) = (-t4./(1-1i.*H1)).*y(3);

in your function definition.

I didn't look further - maybe there are more errors in the setup of your system of equations.

[Syed Mohiuddin]

Thanks a lot, i have made the correction and have got the result. Your every answer is of great help to me. Thank you once again