I've cleaned up your code (see inline comments).
You have to supply the legend with handles that specify what belongs in the legend and you have to store the handle from each iteration.
I've made some guesses in in the legend inputs so if it doesn't work right away, play around with it to get what you want.
hold on ; % MOVE THIS OUTSIDE THE LOOP
col=['r' 'y' 'b' 'g']; % MOVE THIS OUTSIDE THE LOOP
lineHandles = gobjects(1,length(dep(1,:))); % ADD THIS
for i=1:length(dep(1,:))
PP(:,i)=plot(dep(:,i),ext','ok','LineWidth',2)% original data
P1(:,i) = polyfit(dep(:,i),ext,2);%%2 indicates the quadratic
% x_values = min(in_(1,:)): 0.001:max(x_axis(1,:));
Y= polyval(P1(:,i),dep(:,i));
h1 = plot(dep(:,i),Y,'*g','LineWidth',1) % fitted
Y1 = polyval(P1(:,i),dep(:,i));
R1 = corrcoef(Y1,dep(:,i));
%%lsline
ss = min(dep(:,i)) : 0.001: max(dep(:,i));
Y2 = polyval(P1(:,i),ss);
lineHandles(i) = plot(ss,Y2,'color',col(i),'LineWidth',1); % STORE HANDLE
end
% MOVE THIS OUTSIDE THE LOOP
title('Polynomial fitting','FontSize',16,'FontWeight','normal')
xlabel('Depolarisation ratio \delta_{in}','FontSize',12,'FontWeight','normal');
ylabel('Extinction coefficient \alpha{(m^-^1)}','FontSize',12,'FontWeight','normal');
ylim([0.005 0.03])
% SUPPLY HANDLES, DEFINE EACH POLY CURVE LINE
legend([PP(1,i), h1, lineHandles],...
'Actual data','Data fit',compose('PolyCurve %d\mum',1:length(dep(1,:))))
