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