Please help me to run this simple code

%ErrorError using surf (line 71)

X, Y, Z, and C cannot be complex.

Error in Untitled (line 66)

surf(X, Y, Z);

proj

ky
    0.0100

The solution was obtained on a mesh of 754 points.
The maximum residual is  2.042e-06. 
There were 39912 calls to the ODE function. 
There were 274 calls to the BC function. 
The solution was obtained on a mesh of 192 points.
The maximum residual is  8.109e-06. 
There were 6379 calls to the ODE function. 
There were 110 calls to the BC function. 

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

The solution was obtained on a mesh of 2674 points.
The maximum residual is  1.847e-01. 
There were 85426 calls to the ODE function. 
There were 341 calls to the BC function. 

Error using surf (line 71)
X, Y, Z, and C cannot be complex.

Error in solution>proj (line 69)
surf(X, Y, Z);
^^^^^^^^^^^^^^

function sol = proj
    clc; clf; clear;
    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;sigf=0.05*10^-8;
ky=muf/rhof;
disp('ky');
disp((muf/rhof));

%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);
sigt = sigf + (3 * ((ph1 * sig1 + ph2 * sig2) - sigf * (ph1 + ph2)) /(((ph1 * sig1 + ph2 * sig2) / (sigf * (ph1 + ph2))) + 2 - ((ph1 * sig1 + ph2 * sig2) / sigf) + (ph1 + ph2)));
%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)))+(1-ph1)*ph1*((rho1*be1)/(rhof*bef));  
 myLegend1 = {};myLegend2 = {};
 rr = [0.5 1 1.5];
    numn = numel(rr);
    m = linspace(0, 7);
R=3;M=2;Qh=0.1;
    Rd=1;Pr=6.9;
     y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);
   solinit = bvpinit(m,y0);
    
 Z = zeros(numn, length(m));
   
    % Solve the BVP for each Prf
    for i = 1:numn
        n = 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;
    

     
   
    i=i+1;


figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)                                                                                                                                                                                                                                                                               
    hold on
    function dy = projfun(~, y)
      dy= zeros(8,1);
        % alignComments
        f = y(1);
        df = y(2);
         g = y(3);
        dg= y(4);
        h= y(5);
        dh = y(6);
         t = y(7);
       dt=y(8);
               dy(1) = df; 
        dy(2) =(R^((1-n)/2)*(muf/mut)*(rhot/rhof)*(f^2-g^2+R*h*df+(sigt/sigf)*M*f))^(1/n);
         dy(3) = dg;
     dy(4) = (R^((1-n)/2)*(muf/mut)*(rhot/rhof)*(f*g+R*dg+(sigt/sigf)*M*g))^(1/n);
         dy(5) =dh ;
      dy(6) =(R^((-n))*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h))^(1/n);
         dy(7) =dt;
       dy(8)=(1/(1+(4/3)*Rd*(1/(kt/kf))))*(1/R)*(-t-Qh*Pr*(kf/kt)*t+Pr*((vt/(rhof*cpf))*(kf/kt)*(f*t+R*h*dt)));         

         
    end
end

function res= projbc(ya,yb)
res= [ya(1)-0.3;
  ya(3)-1;
  ya(5)-0.1;
    ya(7)-1;
    
    yb(1)-0.3;
        yb(3);
        yb(5);
     yb(7);
       ];
       end

Answer 1

[Torsten]

Your solution for rr(3) becomes complex-valued.

For complex-valued functions, you can plot e.g. their real part, their imaginary part or their absolute value:

surf(X,Y,real(Z))
surf(X,Y,imag(Z))
surf(X,Y,abs(Z))

You try to plot the solution y(1,:). Why do you use a legend that tells you plot y(6,:) ?

Do you really mean

dy(6) =(R^((-n))*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h))^(1/n);

If I look at dy(2) and dy(4), it could be

dy(6) =(R^((-n)*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h)))^(1/n);