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

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