Please help me to run surf plot with bvp4c.

Please help me to run surf plot with bvp4c.The surfce digram consists of (constant Prf [2 :1:6] represents y-axis & vector>> sol.x [0 4] represents x-axis (m = linspace(0,4);) & second solution only sol.y(6,:) represents z-axis).The following is the code for 2D (sol.x [0 4] and only sol.y(6,:)). How to give command for making surf plot in bvp4c.

function sol= proj
clc;clf;clear;
%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 = {};
rr = [4 5 6 7]
for i =1:numel(rr)
    Prf = rr(i);
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;
  p=-0.5; L=(muf/rhof);L1=L^(p);

   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,:)))
   % axis([0 4 0 1])
    grid on,hold on
    myLegend1{i}=['Pr = ',num2str(rr(i))];
    

    
    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) =(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);
        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