wns=[0.1:0.1:1]; z=0.707; yy=[]; t=0:.1:12;
for i=1:length(wns)
   wn=wns(i); 
   y=wn*exp(-z*wn*t).*sin(wn*sqrt(1-z^2)*t)/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,0.5])