Sindar's answer works for R2018b and onward. For 2016a, you can use XLim as the x coordinates of line
clear all;
clc;
close all
x=[33.5 44 55 69 106 143 123 49 -25];
y=[61.8 53 45 73 77 81 156 68 -18];
z=[-7.5 -5 -2.5 -2.5 0 2.5 2.5 5 7.5];
zz=z(1:3:end)
for i=1:3
zz_n(i)=zz(i)
end
ax = axes();
hold on % just one hold on is enough
for i=1:3
plot(x,z)
plot(ax.XLim,zz_n(i)*ones(1,2),':')
plot(y,z)
end
