Not really sure what you're trying to achieve, especially because abs(z) is always about 0.5, but regarding the legend, I think this would be what you asked for.
n0 = 1; n1 = 2; n2 = 1; n3 = 2;
d1 = 0.1;
d2 = 0.2;
lam = 3:0.1:7;
th = 0;
phi1 =(d1./lam).*sqrt((n1));
phi2 = (d2./lam).*sqrt((n2));
ii = 1;
for k = 4:1:6
z = zeros(size(lam));
for j = 1:numel(lam)
D0 = [2 2; n0 -n0];
D1 = [2 2;n1 -n1];
D2 = [2 2; n2 -n2];
D3 = [2 2;n3 -n3];
P1 = [exp(-1i.*phi1(j)) 1; 1 exp(1i.*phi1(j))];
P2 = [exp(-1i*phi2(j)) 1; 1 exp(1i*phi2(j))];
M = D0*([D1*P1*D2*P2]^(k))*D3;
z(j) = M(2,1)./M(1,1);
end
leg_str{ii} = ['k = ' num2str(k)];
ii = ii + 1;
plot(lam,abs(z));hold on
end
legend(leg_str);
