Solution 1:

Firstly, you say that with Laplace transform we can only calculate forced response. I think one can also find the initial condition response with Laplace transform with the formula:

\begin{equation*} L{(\frac{d^n}{dt^n}y(t))}=s^nY(s)-s^{n-1}y(0)-s^{n-2}y'(0)-\cdots -sy^{(n-2)}(0) -y^{(n-1)}(0) \end{equation*}

Secondly, you say that if input contains natural frequencies then Laplace transform fails. I read your answer about Laplace transform, but there you derive general formula by assuming the input frequency is not equal to one of the natural frequencies of the system. In fact, in your formula there you have $H(s_0)$, which is infinity at poles. However, I think, it is that formula that fails, not the Laplace transform. For example, consider your system in this question and input $e^{-2t}$. Then we have:

\begin{equation*} Y(s)=H(s)U(s)=\frac{1}{s^2+3s+2} \frac{1}{s+2}= \frac{1}{(s+2)^2(s+1)} \end{equation*}

Then you can find $y(t)$ as

\begin{equation*} y(t) = L^{-1}(Y(s))= L^{-1}\left(\frac{1}{(s+2)^2(s+1)}\right) = L^{-1}\left(\frac{-1}{(s+2)^2}+\frac{-1}{s+2} + \frac{1}{s+1}\right) \\ = -te^{-2t}-e^{-2t}+e^{-t}. \end{equation*}

Finally, why we prefer Laplace transform? At least in my opinion, Laplace transform is well defined, but operator domain seems little bit more abstract. For example what is the meaning of D in your equations? I think it is a bit abstract. However, Laplace transform is well defined as:

$L(y(t))=\int_0^\infty y(t)e^{-st}dt$

and all the properties follow from the definition. So you straightforwardly understand what is $s$ standing for in the equations.

Any objections or corrections are welcome.