Wrong graph is coming with this code

function main
A=0.5; pr=1; a=1; phi=0.1;rhos=997.1;Cps=4179;ks=0.613;rhof=8933;Cpf=385;kf=401; 
a1=(1-phi)^-2.5*((1-phi)+phi*(rhos/rhof));  Knf=(ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));
a2=((1-phi)+phi*((rhos*Cps)/(rhof*Cpf)))*Knf;
xa = 0;xb = 5;solinit=bvpinit(linspace(xa,xb,101),[0 1 0 a 1 0 0 0]);
sol = bvp4c(@ode,@bc,solinit); x = linspace(xa,xb,101);S = deval(sol,x);
function res = bc(ya,yb)
res = [ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];
end 
function dydx = ode(~,y)    
dydx = [y(2);  y(3);  2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));  
        y(5);  y(6);  2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
        y(8);  A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))];              
end
figure(14)
f1 = @(x,y) S(2);
fsurf(f1)
xlabel '\bfx';ylabel '\bfy';zlabel '\bff^\prime(\eta)'
hold on
end
%% The graph is not obeying Boundary Condition
%% I am trying graphs like Fig. (14)  and  another  is contour (Attached)

[MINATI]

can anyone share idea here

[MINATI]

Dear Star Strider

Please have a try

Answer 1

[Star Strider]

Using fsurf is not appropriate here.

Try this instead:

figure(14)
hs = surf(S);
grid on
xlabel '\bfx';ylabel '\bfy';zlabel '\bff^\prime(\eta)'
% shading('interp')                                               % Optional
% hs.EdgeColor = 'none';                                          % Optional

That produces this plot:

I am not certain what the reference to ‘S(2)’ is suposed to do, beccause ‘S’ is a matrix. This plots the entire matrix.

.

[MINATI]

Dear Star Strider

Thanks for quick reply.

How to draw surf(S(2,:)) with x & y

And also surf(meshgrid(x), meshgrid(phi), S(3,1))

Where phi = linspace(0,0.15,101)

[Star Strider]

My pleasure!

I am completely lost.

With respect to: ‘How to draw surf(S(2,:)) with x & y ’ this works:

hs = surf(x, [2 2], [1;1]*S(2,:))

With respect to: ‘And also surf(meshgrid(x), meshgrid(phi), S(3,1))’ is just not going to work because ‘S(3,1)’ is a point, so plotting it is useless, as welll as the rest of that not making any sense.

With respect to: ‘Where phi = linspace(0,0.15,101)’ note that ‘phi’ is currently a scalar. I have no idea what you are doing, however you will need to incorporate that vector with your code. Every variable that uses ‘phi’ will also become vectors, likely not something bvp4c can handle, since your ODE/BVP function needs to return a column vector.

So the best way of dealing with that may be a loop, looping through the various values of ‘phi’ and storing the results by concatenating the resulting ‘S’ matrices along the third dimension. Then choose the appropriate values (rows, columns, pages) of that 3D matrix for the plots.

When you have that working, post back here with the results. It will likely be possible to plot them.

[Star Strider]

MINATI’s Answer moved here —

Dear Star Strider have a look into the following code:

function main

A = 0.5; pr = 1; %%a = 1; phi = 0.1;

rhos = 997.1;Cps=4179;ks=0.613;rhof=8933;Cpf=385;kf=401;

phiv = linspace (0,0.15,101);

for a = [0 0.2 0.5 1]

for i = 1:length(phiv)

phi = phiv(i);

a1=(1-phi)^2.5*((1-phi)+phi*(rhos/rhof)); Knf=(ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));

a2 =((1-phi)+phi*((rhos*Cps)/(rhof*Cpf)))*Knf;

BCres = @(ya,yb)[ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];

fODE = @(x,y)[y(2); y(3); 2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));

y(5); y(6); 2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));

y(8); A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))];

xa = 0;xb = 5;solinit=bvpinit(linspace(xa,xb,101),[0 1 0 a 1 0 0 0]);

sol = bvp4c(fODE,BCres,solinit);x = linspace(xa,xb,101);S = deval(sol,x);

Sk(i,1) = S(3,1)/((1-phi)^2.5); Nu(i,1) = - Knf * S(5,1);

end

figure(10)

plot(phiv,Sk,'-','Linewidth',1.5)

% % title '(b)'

xlabel('\bf\phi');

ylabel('\bfSk_{r}')

hold on

end

end

%%%%%%%%I have incorporated a loop for phi variation where plotting of S(3,1) is OK but

surf(S(3,:)) is not coming

Can it possible to modify for surf of Sk (OR In x-axis, phi variation In Y-axis for a = linspace(0,1,101))

Please do something

[Star Strider]

Create a second loop for ‘a’:

function main
A = 0.5; pr = 1; %%a = 1; phi = 0.1;
rhos = 997.1;Cps=4179;ks=0.613;rhof=8933;Cpf=385;kf=401; 
phiv = linspace (0,0.15,101);
av = [0 0.2 0.5 1];
for k = 1:numel(av)
    a = av(k);
 for i = 1:length(phiv)
     phi = phiv(i);
a1=(1-phi)^2.5*((1-phi)+phi*(rhos/rhof));  Knf=(ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));
a2 =((1-phi)+phi*((rhos*Cps)/(rhof*Cpf)))*Knf;
BCres = @(ya,yb)[ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];  
fODE = @(x,y)[y(2);  y(3);  2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));  
        y(5);  y(6);  2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
        y(8);  A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))]; 
xa = 0;xb = 5;solinit=bvpinit(linspace(xa,xb,101),[0 1 0 a 1 0 0 0]);
sol = bvp4c(fODE,BCres,solinit);x = linspace(xa,xb,101);S = deval(sol,x);
Sk(i,k) = S(3,1)/((1-phi)^2.5);   Nu(i,1) =  - Knf * S(5,1);  
 end
end

figure(10)
surf(av, phiv, Sk)
% % title '(b)'
xlabel('a')
ylabel('\bf\phi'); 
zlabel('\bfSk_{r}')
view(40,40)
end

Also, put the surf call after both loops.

This then creates a matrix for ‘Sk’ that you can use with surf.

.

[MINATI]

Thanks Star Strider

You are awsome.

Any way the first problem remain unanswered (Attached pdf)

surf(x, y, S(2,:)) and contour

pdf exist means code can be created but I can't.

But I have confidence on You in this matter. If you agree take your time to create.

A great thanks!

[Star Strider]

Thank you!

Your ccode does not produce the plot in the .PDF file. I have no idea what you are doing.

You can use essentially the same code to create the contour plot as the surf plot. Just call contour instead of surf, and remove the zlabel and view calls.

As always, my pleasure!