The loop is not necessary.
This required a few corrections (a typographical error and to vectorise the exponentiation), and now appears to work:
Es=200000;
Esh=8500;
fy=500;
fsu=750;
epssh=0.009;
epssu=0.075;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 10, 500)*1E-3;
sigmaSteel=@(epscm) Es*epscm .* (epscm<=epsy) + fy .* (epscm>epsy & epscm<=epssh) + fsu+(fy-fsu)*abs((epssu-epscm)./(epssu-epssh)).^(1/P) .* (epscm>epssh & epscm<=epssu) + 0 .* (epscm>epssu);
figure
plot(epscmv, sigmaSteel(epscmv))
grid on
It does not exactly reproduce the plot in the 1.jpg attachment, although it appears to be reasonably close.
Best Answer