How to format the figure from 2-D to surface figure like as the attached figure. The third axis represnt different values for Bi "Take any value for Bi"

Question

0 votes

function sol= proj
clc;clf;clear;
global lamda
%Relation of base fluid
rhof=1;
kf=0.613*10^5;
cpf=4179*10^4;
muf=10^-3*10;
sigf=0.05*10^-8;
alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;
rho1=5180*10^-3;
cp1=670*10^4;
k1=9.7*10^5;
sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=401*10^5;
sig2=5.96*10^-1;
%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};myLegend2 = {};
rr = [1 10 20]
for i =1:numel(rr)
    Re = rr(i);
    Bi=1;
    
    Prf=2;
p=-0.5; L=(muf/rhof);L1=L^(p);Lt=L1^+1;
Nr=0.1;
    gamma=pi/3;
a=1;b=0.1;v=1;u=1;
   M=3;
  Nt=1;Nb=1; sc=0.6;s1=1;s2=1;
  
   Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
   Lf=rhof*kf;
   y0 = [1,0,1,0,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,4);
   solinit = bvpinit(m,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
      disp((sol.y(1,20)))
    figure(1)
    plot(sol.x,(sol.y(6,:)))
    grid on,hold on
    myLegend1{i}=['Pr = ',num2str(rr(i))];
xlabel('eta');
    ylabel('(thetas-thetaf)/thetas');
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on
    function dy= projfun(~,y)
        dy= zeros(9,1);
%         alignComments
        E  = y(1);
        dE = y(2);
        F  = y(3);
        dF= y(4);
        W = y(5);
        t = y(6);
        dt  = y(7);
        phi = y(8);
        dphi  = y(9);
       
    dy(1) = dE;        
          dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
         dy(3) = dF;
          dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));
         dy(5) = -(1/L1)*(u*a*E+b*v*F);
        dy(6) = dt;
dy(7) =(Bi*Prf*(rhof/muf)*Re)*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
        dy(8)= dphi;
dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
    
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5)-0;
    ya(6)-1;
    ya(8)-1;
    yb(1);
    yb(3);
    yb(6);
    yb(8)];
end