please help me to add legend in run figure as attached photo

function sol= proj
clc;clf;clear;
myLegend1 = {};
rr = [0 0.5 1 1.5 2]
numR = numel(rr);
    m = linspace(0,30);
    alfa=0.6;
     L1=2;L2=-0.1;
    
    s=0.01;h=0.02;
    
  k0=386;  ce=3.831*10^2;   mu=38.6*10^9;alfat=1.78*10^-5;  rho=89.54*10^2;   lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
        c0=sqrt((lamda+2*mu)/(rho));  Betanot=(3*lamda+2*mu)*alfat;  a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
       a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
    a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);

Z = zeros(numR, length(m)); y0 = [0,0,0,1];
options =bvpset('stats','on','RelTol',1e-5);
for i =1:numR
    t= rr(i);
    
    if i==1
      solinit = bvpinit(m, y0);
    else
      guess = @(x)interp1(sol.x,(sol.y).',x);
      solinit = bvpinit(sol.x,guess);
    end
    sol = bvp4c(@projfun,@projbc, solinit, options);
    Z(i,:) = interp1(sol.x,sol.y(1,:),m);
    i=i+1;
end
[X, Y] = meshgrid(m, rr);
    shading interp;
    legend('hh')
   figure(1)
surf(X, Y, Z);
grid on,hold on
    myLegend1{i}=['alfa= ',num2str(rr(i))];
    figure(1)
legend(myLegend1)
hold on
xlabel('x');
ylabel('t');
zlabel('T');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
      
%       N1=zeros(1,1); 
%             N1=(s*L1+a4*h*L2)*sol.y(2,:)-alfa*(s*L1+a4*h*L2)*sol.y(1,:).*sol.y(2,:)-a3*sol.y(1,:)+a3*alfa*sol.y(1,:).^2;
%            figure(1)
%             plot(sol.x,N1)
%              title('sigma xx')
%                 grid on,hold on
%                 myLegend1{i}=['alfa= ',num2str(rr(i))];
% figure(2)
%    plot(sol.x,(sol.y(1,:)))
%    title('Temperature')
%    grid on,hold on
%                    myLegend2{i}=['alfa= ',num2str(rr(i))];
% 
%     figure(3)
%     title('Temperature gradient')
%      plot(sol.x,(sol.y(2,:)))
%     grid on,hold on
%     myLegend3{i}=['alfa= ',num2str(rr(i))];
%     
%        figure(4)
%     plot(sol.x,L1*(sol.y(1,:)))
%      title('U')
%     grid on,hold on
%     myLegend4{i}=['alfa= ',num2str(rr(i))];
%      figure(5)
%      title('V')
%      plot(sol.x,L2*(sol.y(1,:)))
%     grid on,hold on
%     myLegend5{i}=['alfa= ',num2str(rr(i))];
% 
%      N2=zeros(1,1); 
%             N2=a1*(h*L1+s*L2)*(sol.y(2,:)-alfa*sol.y(1,:).*sol.y(2,:));
%            figure(6)
%             plot(sol.x,N2)
%              title('sigma xy')
%                 grid on,hold on
%                 myLegend6{i}=['alfa= ',num2str(rr(i))];
%     

    function dy= projfun(x,y)
        dy= zeros(4,1);
        E  = y(1); 
        dE = y(2);
        ddE=y(3);
        dddE=y(4);
   dy(1) = dE;  
   dy(2) =ddE;
   dy(3)=dddE;
   dy(4)=(1/((L1^2*(s^2+a1*h^2-(x+t)^2)+L2^2*(a2*s*h))+(-(alfa*(s^2+a1*h^2-(x+t)^2))*L1^2-(alfa*(a2*s*h))*L2^2)*E))*(-(-(2*(x+t))*L1-a3*s)*dddE+((alfa*2*(x+t))*L1+2*alfa*a3*s)*E*dddE-3*(alfa*2*(x+t)*L1+2*alfa*a3*s)*dE*ddE);
     

    end
    end

function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(2);
    yb(1);
    yb(2);
  ];
end

Answer 1

[dpb]

function sol= proj
clc;clf;clear;
myLegend1 = {};
rr = [0 0.5 1 1.5 2];
numR = numel(rr);
m = linspace(0,30);
alfa=0.6;
L1=2;L2=-0.1;
s=0.01;h=0.02;
k0=386;  ce=3.831*10^2;   mu=38.6*10^9;alfat=1.78*10^-5;  rho=89.54*10^2;   lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
c0=sqrt((lamda+2*mu)/(rho));  Betanot=(3*lamda+2*mu)*alfat;  a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
Z = zeros(numR, length(m)); y0 = [0,0,0,1];
options =bvpset('stats','on','RelTol',1e-5);
for i =1:numR
    t= rr(i);
    if i==1
      solinit = bvpinit(m, y0);
    else
      guess = @(x)interp1(sol.x,(sol.y).',x);
      solinit = bvpinit(sol.x,guess);
    end
    sol = bvp4c(@projfun,@projbc, solinit, options);
    Z(i,:) = interp1(sol.x,sol.y(1,:),m);
    i=i+1;
end
[X, Y] = meshgrid(m, rr);
shading interp;
%legend('hh')
figure(1)
surf(X, Y, Z);
grid on,hold on
%myLegend1{i}=['alfa= ',num2str(rr(i))];
myLegend1=sprintf('alfa= %g',alfa);
legend(myLegend1)
hold on
xlabel('x');
ylabel('t');
zlabel('T');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
      
%       N1=zeros(1,1); 
%             N1=(s*L1+a4*h*L2)*sol.y(2,:)-alfa*(s*L1+a4*h*L2)*sol.y(1,:).*sol.y(2,:)-a3*sol.y(1,:)+a3*alfa*sol.y(1,:).^2;
%            figure(1)
%             plot(sol.x,N1)
%              title('sigma xx')
%                 grid on,hold on
%                 myLegend1{i}=['alfa= ',num2str(rr(i))];
% figure(2)
%    plot(sol.x,(sol.y(1,:)))
%    title('Temperature')
%    grid on,hold on
%                    myLegend2{i}=['alfa= ',num2str(rr(i))];
% 
%     figure(3)
%     title('Temperature gradient')
%      plot(sol.x,(sol.y(2,:)))
%     grid on,hold on
%     myLegend3{i}=['alfa= ',num2str(rr(i))];
%     
%        figure(4)
%     plot(sol.x,L1*(sol.y(1,:)))
%      title('U')
%     grid on,hold on
%     myLegend4{i}=['alfa= ',num2str(rr(i))];
%      figure(5)
%      title('V')
%      plot(sol.x,L2*(sol.y(1,:)))
%     grid on,hold on
%     myLegend5{i}=['alfa= ',num2str(rr(i))];
% 
%      N2=zeros(1,1); 
%             N2=a1*(h*L1+s*L2)*(sol.y(2,:)-alfa*sol.y(1,:).*sol.y(2,:));
%            figure(6)
%             plot(sol.x,N2)
%              title('sigma xy')
%                 grid on,hold on
%                 myLegend6{i}=['alfa= ',num2str(rr(i))];
%     

    function dy= projfun(x,y)
        dy= zeros(4,1);
        E  = y(1); 
        dE = y(2);
        ddE=y(3);
        dddE=y(4);
   dy(1) = dE;  
   dy(2) =ddE;
   dy(3)=dddE;
   dy(4)=(1/((L1^2*(s^2+a1*h^2-(x+t)^2)+L2^2*(a2*s*h))+(-(alfa*(s^2+a1*h^2-(x+t)^2))*L1^2-(alfa*(a2*s*h))*L2^2)*E))*(-(-(2*(x+t))*L1-a3*s)*dddE+((alfa*2*(x+t))*L1+2*alfa*a3*s)*E*dddE-3*(alfa*2*(x+t)*L1+2*alfa*a3*s)*dE*ddE);
     

    end
    end

function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(2);
    yb(1);
    yb(2);
  ];
end

proj;

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

The solution was obtained on a mesh of 3592 points.
The maximum residual is  2.858e-01. 
There were 150449 calls to the ODE function. 
There were 298 calls to the BC function. 
The solution was obtained on a mesh of 1722 points.
The maximum residual is  8.459e-07. 
There were 27575 calls to the ODE function. 
There were 19 calls to the BC function. 
The solution was obtained on a mesh of 1722 points.
The maximum residual is  3.242e-08. 
There were 34435 calls to the ODE function. 
There were 32 calls to the BC function. 
The solution was obtained on a mesh of 1722 points.
The maximum residual is  4.936e-09. 
There were 34434 calls to the ODE function. 
There were 32 calls to the BC function. 
The solution was obtained on a mesh of 1722 points.
The maximum residual is  2.011e-09. 
There were 27546 calls to the ODE function. 
There were 19 calls to the BC function. 

●

There was a superfluous call to legend

legend('hh')

I commented out; I presume that may have just been an attempt to see if a legend would show up at all that wasn't then removed.

The real problem is that line creating the legend text is outside (after) the loop on i over the values in rr so the reference to rr(i) is out of bounds in the original code. It wouldn't have provided the desired value of the variable alfa, anyway.

myLegend1{i}=['alfa= ',num2str(rr(i))];

What you were looking for is

myLegend1=sprintf('alfa= %g',alfa);
legend(myLegend1)

instead.