The loop is not necessary.
This required a few corrections (a typographical error and to vectorise the exponentiation), and now appears to work:
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);
It does not exactly reproduce the plot in the 1.jpg attachment, although it appears to be reasonably close.