wn=1;
more on
sigma=[0,0.2,0.4,0.6,0.9,1.2,1.5];
num=wn*wn
t=linspace(0,20,200)'
for j=1:7;
    den=conv([1,0],[1,2*wn*sigma(j)]);
    s1=tf(num,den);
    sys=feedback(s1,1)
    y(:,j)=step(sys,t);
    %y=step(sys,t);
    %hold on
    % plot(t,y);
end
 

   plot(t,y(:,1:7));