Please help me to run this code

function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;

%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;

%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
%for i =1:numel(rr)
rr = [0 1 2];
numGr = numel(rr);
m = linspace(0,1);
a=-0.001;b=0.0001;p=-0.15/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=0.5;
gamma=pi/4;
prf=6.9;Rd=0.5;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
disp('coe');disp((mut/muf)*(rhof/rhot));
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
%solinit = bvpinit(m,y0);
%  sol= bvp4c(@projfun,@projbc,solinit,options);
Z = zeros(numGr, length(m));
for i = 1:numGr
    Gr= rr(i);
    solinit = bvpinit(m, y0);
    sol = bvp4c(@projfun, @projbc, solinit, options);
    Z(i, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
figure;
surf(X, Y, Z);
xlabel('x');
ylabel('Prf');
zlabel('Solution y(6,:)');
title('Surface Plot of Solution');
grid on;

    function dy= projfun(~,y)
    	dy= zeros(8,1);
    	%         alignComments
    	E  = y(1);
    	dE = y(2);
    	F  = y(3);
    	dF= y(4);
    	w = y(5);
    	dw=y(6);
    	t = y(7);
    	dt  = y(8);
    	dy(1) = dE;
    	dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
    	dy(3) = dF;
    	dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
    	dy(5) =-(a*F+b*E);
    	dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+Gr*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
    	dy(7) = dt;
    	dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
    	
    end
end

function res= projbc(ya,yb)
    
    res= [ya(1)+1;
    	ya(3)-1;
    	ya(5);
    	ya(6);
    	ya(7)-1-(1/0.9)*ya(8);
    	yb(1)-0.01;
    	yb(3);
    	yb(7);
    	% yb(7);
    	];
end

Answer 1

[Torsten]

Note that you plot y(1,:), not y(6,:).

If I were you, I'd simply use "plot" to compare the three curves depending on Gr.

proj()

ky
    0.0100

coe
    0.8935

The solution was obtained on a mesh of 100 points.
The maximum residual is  2.908e-07. 
There were 2614 calls to the ODE function. 
There were 54 calls to the BC function. 
The solution was obtained on a mesh of 100 points.
The maximum residual is  3.207e-07. 
There were 2895 calls to the ODE function. 
There were 55 calls to the BC function. 
The solution was obtained on a mesh of 100 points.
The maximum residual is  2.168e-06. 
There were 2895 calls to the ODE function. 
There were 55 calls to the BC function. 

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909 0.1010 0.1111 0.1212 0.1313 0.1414 0.1515 0.1616 0.1717 0.1818 0.1919 0.2020 0.2121 … ] (1×100 double)
         y: [8×100 double]
        yp: [8×100 double]
     stats: [1×1 struct]

function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;

%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;

%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
%for i =1:numel(rr)
rr = [0 1 2];
numGr = numel(rr);
m = linspace(0,1);
a=-0.001;b=0.0001;p=-0.15/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=0.5;
gamma=pi/4;
prf=6.9;Rd=0.5;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
disp('coe');disp((mut/muf)*(rhof/rhot));
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
%solinit = bvpinit(m,y0);
%  sol= bvp4c(@projfun,@projbc,solinit,options);
Z = zeros(numGr, length(m));
for i = 1:numGr
    Gr= rr(i);
    solinit = bvpinit(m, y0);
    sol = bvp4c(@projfun, @projbc, solinit, options);
    Z(i, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
surf(X, Y, real(Z));
xlabel('x');
ylabel('Prf');
zlabel('Solution y(6,:)');
title('Surface Plot of Solution');
grid on

    function dy= projfun(~,y)
    	dy= zeros(8,1);
    	%         alignComments
    	E  = y(1);
    	dE = y(2);
    	F  = y(3);
    	dF= y(4);
    	w = y(5);
    	dw=y(6);
    	t = y(7);
    	dt  = y(8);
    	dy(1) = dE;
    	dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
    	dy(3) = dF;
    	dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
    	dy(5) =-(a*F+b*E);
    	dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+Gr*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
    	dy(7) = dt;
    	dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
    	
    end
end

function res= projbc(ya,yb)
    
    res= [ya(1)+1;
    	ya(3)-1;
    	ya(5);
    	ya(6);
    	ya(7)-1-(1/0.9)*ya(8);
    	yb(1)-0.01;
    	yb(3);
    	yb(7);
    	% yb(7);
    	];
end

[Tarek]

Thanks for your positive feedback.

Finally, please i want to draw above figure as the position of the attached figure

[Tarek]

Here you are professor Toreston

[Torsten]

Finally, please i want to draw above figure as the position of the attached figure

I don't understand what you mean.

[Tarek]

I want to make the 3D figure from the above code inclined as you like the attached photo

[Torsten]

The picture should pop up automatically when you run the code. But looking at your recent question in the forum, you seem to have problems with the "surf" function.