wns=[0.1:0.1:1]; z=0.55; yy=[]; t=0:.1:12;
for i=1:length(wns)
   wn=wns(i); %z in (0,1) 
   wd=wn*sqrt(1-z^2); th=atan(sqrt(1-z^2)/z);
   y=1-exp(-z*wn*t).*sin(wd*t+th)/sqrt(1-z^2);
   yy=[yy;y];
end
plot(t,yy)

mesh(wns,t,yy'); 
set(gca,'XDir','reverse'); axis([0,1,0,12,0,1.2])