Could anyone help me with my matlab bvp4c program?

Hello Dear Sir/ Medam I have two couppled equations (temparature and Concentration) with Boundary conditions. I have solved this equations and solved analytically and got the graph. Now i have used BVP4C for the same equations and got graph is not maching with my earlier graph. Both program and graph i am sharing hear, if u help me to find where i got wrong. It will be helpful to me. Thank you.

Answer 1

[Torsten]

The order of your functions in the dydy function handle is wrong.

It should work now.

%I am a Researcher, my problem is  i try to solve  couppled ODE s with BVP4C MATLAB code the graph obtained after running the MATLAB program, which differs from the actual graph obtained using Mathematica software. Here i am sharing the Equations and Boundary Conditions.
% temp Eqn: Theta^''=-N1*Pr*G6* (dTheta/dy)+S1*Pr*G7*Theta;
% Concen Eqn: Phi^''= -N1*Sc*G8*(dPhi/dy)+H1*Sc*G8*Phi;
% Boundary conditions: dTheta/dy=G7*B1*( Theta-1),dPhi/dy=G8*F1*( Phi-1)  at y=0.
%                     Theta,Phi = 0 at y tends to infinity.
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.15;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
dydx = @(x,y)[y(2);
             -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T1 = bvp4c(dydx, BC, solint);
hold on

plot(T1.x,([T1.y(1,:)]),'b','linewidth',1.5)

B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.25;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T2 = bvp4c(dydx, BC, solint);
hold on

plot(T2.x,([T2.y(1,:)]),'r','linewidth',1.5)


B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.35;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T3 = bvp4c(dydx, BC, solint);
hold on

plot(T3.x,([T3.y(1,:)]),'g','linewidth',1.5)

[sunitha kolasani]

Thank you very much.

[Ramanuja]

Thank you very much. send the question sir

[sunitha kolasani]

How can i smooth above curve . Is there any suggestions?

[Torsten]

Use more points for plotting in x-direction:

%I am a Researcher, my problem is  i try to solve  couppled ODE s with BVP4C MATLAB code the graph obtained after running the MATLAB program, which differs from the actual graph obtained using Mathematica software. Here i am sharing the Equations and Boundary Conditions.
% temp Eqn: Theta^''=-N1*Pr*G6* (dTheta/dy)+S1*Pr*G7*Theta;
% Concen Eqn: Phi^''= -N1*Sc*G8*(dPhi/dy)+H1*Sc*G8*Phi;
% Boundary conditions: dTheta/dy=G7*B1*( Theta-1),dPhi/dy=G8*F1*( Phi-1)  at y=0.
%                     Theta,Phi = 0 at y tends to infinity.
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.15;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
dydx = @(x,y)[y(2);
             -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T1 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T1);
plot(xplot,y1plot(1,:),'b','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.25;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T2 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T2);
plot(xplot,y1plot(1,:),'r','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.35;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
%              y(4);
%              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
%              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
              -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
              y(4);
              -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
              yb(1);
              ya(4)- F1.*G8*(ya(3)-1);
              yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T3 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T3);
plot(xplot,y1plot(1,:),'g','linewidth',1.5)

%plot(T3.x,([T3.y(1,:)]),'g','linewidth',1.5)

[sunitha kolasani]

Thank you very much its work for me.