Please help me....., the masseage error at the

%Error using surf

%Too many input arguments.

%Error in proj (line 70)

%surf(X, Y, real(Z));

%code
function sol= proj
clc;clf;clear;
global n;
%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 = {};


a=0.009;b=0.45;c=12;p=0.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Gr=0.01*bef^2; 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];
rr = [0.5 1 1.5];
m = linspace(0,1);
numn = numel(rr);
options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numn, length(m));

for K = 1:numn
    n= rr(K);
    solinit = bvpinit(m, y0);
    sol = bvp4c(@projfun, @projbc, solinit, options);
    Z(K, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, real(Z));
xlabel('eta');
ylabel('Gr');
zlabel('z');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;


    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)*(c*n*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*c*n*sin(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(3) = dF;
        dy(4) = (((rhot/mut)*(c*n*(muf/rhof)*w*dF)))^(1/(n));
        dy(5) =-(a*F+b*E);
        dy(6) = (((rhot/mut)*((c*n)^(1-n)*(muf/rhof)^(-0.5)*w*(-(a*F+b*E))+Gr*(c*n)^(1-n)*(muf/rhof)^(-0.5)*cos(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(7) = dt;
        dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*(muf/rhof)*((vt/(rhof*cpf))*w*dt-(mut/muf)*Ec*(1/(c*n))*(dw)^2*(muf/rhof)^(1/(2*n))) ;   
        
    end
end
function res= projbc(ya,yb)

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

[Stephen23]

My guess is that you have shadowed SURF with some other function. What is the output of this command?:

which surf -all

/MATLAB/toolbox/matlab/graphics/graphics/graph3d/surf.m
surf is a built-in method                                % tall method

In any case:

proj()

ky
    0.0100

coe
    0.8935

The solution was obtained on a mesh of 145 points.
The maximum residual is  7.907e-06. 
There were 15153 calls to the ODE function. 
There were 163 calls to the BC function. 
The solution was obtained on a mesh of 100 points.
The maximum residual is  6.486e-06. 
There were 3694 calls to the ODE function. 
There were 79 calls to the BC function. 

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 1006 points and the solution are available in the output argument.
The maximum residual is 0.00256409, while requested accuracy is 1e-05.

The solution was obtained on a mesh of 1513 points.
The maximum residual is  2.564e-03. 
There were 308348 calls to the ODE function. 
There were 1122 calls to the BC function. 

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 1.2830e-07 2.5659e-07 3.8489e-07 5.1318e-07 1.9886e-06 3.4640e-06 5.1960e-06 6.9280e-06 8.9487e-06 1.0969e-05 1.3567e-05 1.6165e-05 1.8282e-05 … ] (1×1006 double)
         y: [8×1006 double]
        yp: [8×1006 double]
     stats: [1×1 struct]

function sol = proj()
global n;
%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 = {};
a=0.009;b=0.45;c=12;p=0.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Gr=0.01*bef^2; 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];
rr = [0.5 1 1.5];
m = linspace(0,1);
numn = numel(rr);
options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numn, length(m));

for K = 1:numn
    n= rr(K);
    solinit = bvpinit(m, y0);
    sol = bvp4c(@projfun, @projbc, solinit, options);
    Z(K, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, real(Z));
xlabel('eta');
ylabel('Gr');
zlabel('z');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
    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)*(c*n*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*c*n*sin(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(3) = dF;
        dy(4) = (((rhot/mut)*(c*n*(muf/rhof)*w*dF)))^(1/(n));
        dy(5) =-(a*F+b*E);
        dy(6) = (((rhot/mut)*((c*n)^(1-n)*(muf/rhof)^(-0.5)*w*(-(a*F+b*E))+Gr*(c*n)^(1-n)*(muf/rhof)^(-0.5)*cos(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(7) = dt;
        dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*(muf/rhof)*((vt/(rhof*cpf))*w*dt-(mut/muf)*Ec*(1/(c*n))*(dw)^2*(muf/rhof)^(1/(2*n))) ;   
        
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5);
    ya(6);
    ya(7)+1-(1/0.01)*ya(8);
    yb(1)-0.01;
    yb(3);
    yb(7);
    % yb(7);
    ];
end

[Tarek]

Hello proffesor Stephen23.

Please help me...

Iam coppied the above code, but this message error showed

%Error using surf

Too many input arguments.

Error in Untitled (line 64)

surf(X, Y, real(Z));%

% code
which surf -all
function sol = proj()
%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 = {};
a=0.009;b=0.45;c=12;p=0.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Gr=0.01*bef^2; 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];
rr = [0.5 1 1.5];
m = linspace(0,1);
numn = numel(rr);
options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numn, length(m));

for K = 1:numn
    n= rr(K);
    solinit = bvpinit(m, y0);
    sol = bvp4c(@projfun, @projbc, solinit, options);
    Z(K, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, real(Z));
xlabel('eta');
ylabel('Gr');
zlabel('z');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
    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)*(c*n*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*c*n*sin(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(3) = dF;
        dy(4) = (((rhot/mut)*(c*n*(muf/rhof)*w*dF)))^(1/(n));
        dy(5) =-(a*F+b*E);
        dy(6) = (((rhot/mut)*((c*n)^(1-n)*(muf/rhof)^(-0.5)*w*(-(a*F+b*E))+Gr*(c*n)^(1-n)*(muf/rhof)^(-0.5)*cos(gamma)*(vb/(rhof*bef))*t)))^(1/(n));
        dy(7) = dt;
        dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*(muf/rhof)*((vt/(rhof*cpf))*w*dt-(mut/muf)*Ec*(1/(c*n))*(dw)^2*(muf/rhof)^(1/(2*n))) ;   
        
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5);
    ya(6);
    ya(7)+1-(1/0.01)*ya(8);
    yb(1)-0.01;
    yb(3);
    yb(7);
    % yb(7);
    ];
end

[Stephen23]

Run this line:

which surf -all

show us what it displays.

[Image Analyst]

Do that on the command window, not in your code. You put it in your code, which is OK, but you have to remember to go to the command window to see how many surf functions you have. Most likely you have written your own surf function which is being called instead of the built-in one.