wn = 2.5;
xi = 0.4;
t = [0:0.05:6];
t1 = acos(xi)*ones(1,length(t));
a1 = (1/sqrt(1-xi^2));
hl = 1 - a1 * exp( -xi * wn * t) .* sin( wn * sqrt(1-xi^2) * t + t1);
bu = a1 * exp( -xi * wn * t) + 1;
bl = 2 - bu;
plot(t,hl,'k-',t,bu,'-.',t,bl,':',t,ones(size(t)),'--');
xlabel('\omega_nt'),ylabel('h(t)');grid on;