wn=1;
num=[0 0 wn^2];

zeta=0:0.25:1.5;
for i=1:1:length(zeta)
    e=zeta(i);
    den=[1 2*e*wn wn^2]
    sys=tf(num,den);
    step(sys,0:0.001:20);
    hold on
    pause
end

hold off