How do I use 'legend' to show for which parameter value is for which graph?

The equations I've used are :

 f"=g(g^2+gamma^2)/(g^2+lambda*gamma^2) ------ (1)
g'= 1/3f'^2-2/3(f*f")+ Mn*f' ------------------------(2)
t"+ 4*Rd*t"/3+ 2*Pr*f*t'/3+ Nb*t'*p'+Nt*(t')^2= 0------(3)
p"+(2*Lew*f*p')/3+ Nt*t"/Nb= 0 ------------------------(4)

Where 'lambda', 'gamma', 'Rd', 'Pr', 'Lew', 'Nb', 'Nt', 'Mn' are some parameters
where the boundary conditions are : f=0, f'= delta, t=1, p=1 at y=0 f'-->0, t-->0, p-->0 as y----> infinity

I have used bvp4c to solve these above equations, I have plotted the graphs varying the values of 'lambda' by a 'for' loop. Please, could someone show how do I use the 'legend' feature to show which graph is for which parameter value?

The 'legend' command that I have used doesn't exactly show which graph is for which parameter using color lines.

 function sol= proj
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt delta Mn
delta=-1;
gama=1;
Pr=2;
Lew=2; 
Nb=1;
Nt=2;
Mn= 1;
pp=[0.5:0.5:1.5];
Rd=2;
figure(1)
plot(4,1);hold on
legend('lambda')
options=bvpset('stats','on','RelTol',1e-9);
solinit= bvpinit([0:0.01:10],[1,0,0,0,0,0,0]);
for lambda= pp
    sol= bvp4c(@projfun,@projbc,solinit,options);
solinit= sol;
plot(2,1);plot(sol.x,sol.y(6,:));hold on
legend(['lambda=',num2str(pp)]);
end
end
function f= projfun(x,y)
global lambda gama Pr Rd Lew Nb Nt Mn
f= [y(2)
    y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2)
    y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))+Mn*y(2)
    y(5)
    -((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/(1+4*Rd/3)
    y(7)
    -(2*Lew*y(1)*y(7))+ Nt*((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/Nb*(1+4*Rd/3)];
end
function res= projbc(ya,yb)
global delta
res= [ya(1); ya(2)-delta; ya(4)-1.0; ya(6)-1.0; yb(2); yb(4); yb(6)];
end