X=[]; Y=[]; w=logspace(-1,1,100);
for i=1:length(zeta)
   dd=wn^4+w.^4+2*(2*zeta(i)^2-1)*wn^2*w.^2;
   x=wn^2*(wn^2-w.^2)./dd;
   y=-2*zeta(i)*wn^3*w./dd;
   X=[X, x']; Y=[Y, y'];
end
plot(X,Y), axis('square')
xlabel('Real Axis'); ylabel('Imag Axis')