syms t s;

sigma = heaviside(t);

SIGMA(s) = laplace(sigma,s);

syms s D w_0 w_e;

assume(w_0>0);
assumeAlso(0<D<1);

G(s) = 1/((1 / w_0^2)*s^2+(2*D/w_0)*s+1);

H(s) = SIGMA * G;

h(t) = ilaplace(H, t);

w_e = w_0*(1 - D^2)^(1/2);

disp(simplify(h,"Steps",100));