wn=1; zetas=[0:0.1:1,2,3,5]; yy=[]; t=0:.1:12;
for i=1:length(zetas)
   z=zetas(i);
   if z==0, y=wn*sin(wn*t);
   elseif (z>0 & z<1), 
      y=wn*exp(-z*wn*t).*sin(wn*sqrt(1-z^2)*t)/sqrt(1-z^2);
   elseif z==1, y=wn^2*t.*exp(-wn*t);
   elseif z>1, 
      lam1=-z-sqrt(z^2-1); lam2=-z+sqrt(z^2-1); 
      y=0.5*wn*(exp(lam2*t)-exp(lam1*t))/sqrt(z^2-1);
   end
   yy=[yy;y];
end
plot(t,yy)