How to obtain the solution of a differential equation using a convolution integral? I need to express the solution of this initial value problem about vibration below using Convolution Integral;

$$my''+cy'+ky=f(t) \quad y(0)=0,\quad y'(0)=0$$

But don't have any idea where do i use the Convolution Integral. So how do I do it?
I tried to take laplace transform of both sides.
$$
(ms^2+cs+k)Y(s)=L(f(t))\quad (assuming \quad L(y(t))=Y(s))
$$
I presumed;
$$
f(t)=\int_0^t g(t-T)h(T) \,dT \quad L(g(t))=G(s) \quad L(h(t))=H(s)
$$
So the Convolution theorem gives me the laplace transform of the right side;
$$
L(f(t))=G(s)H(s)
$$
and putting it into the equation;
$$
(ms^2+cs+k)Y(s)=G(s)H(s)
$$
$$
Y(s)=\frac{G(s)H(s)}{ms^2+cs+k}
$$
$$
y(t)=L^{-1}(\frac{G(s)H(s)}{ms^2+cs+k})
$$
I don't know if this solution is enough or correct.
 A: Method by Laplace transform (zero initial conditions):
$$my''+cy'+ky=f(t) \implies [ms^2+cs+k]y(s)=f(s)$$
$$y(s)=\frac{1}{m}\frac{1}{s^2+\frac{c}{m}s+\frac{k}{m}}f(s)$$
The inverse Laplace transform of the last expression is given as
$$y(t)=\frac{1}{m}g(t)*f(t)=\frac{1}{m}\int_{0}^{t}g(t-u)f(t)du$$
In which $g(t)=\mathcal{L}^{-1}\left[\frac{1}{s^2+\frac{c}{m}s+\frac{k}{m}} \right]=\frac{1}{a-b}\left[\exp(at)-\exp(bt) \right]$ and $a$ and $b$ are the roots of $ s^2+\frac{c}{m}s+\frac{k}{m}$.
So we obtain:
$$y(t)=\frac{1}{m}\int_{0}^{t}g(t-u)f(u)du=\frac{1}{m(a-b)}\int_{0}^{t}\left[\exp(at-au)-\exp(bt-bu) \right]f(u)du$$
Solution in state space form:
First cast the system into first-order form:
$$\dot{x}_1 = x_2$$
$$\dot{x}_2 = -x_1-2x_2+f(t)$$
Note the system matrix is given by:
$$\boldsymbol{A}=\begin{bmatrix}0 & 1 \\ -1 &-2\end{bmatrix}$$ 
The input matrix/vector is given by:
$$\boldsymbol{B} =\begin{bmatrix}0\\1\end{bmatrix}$$
Using these expressions we can rewrite the system as
$$\boldsymbol{\dot{x}}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{B}f(t)$$
The solution of such a system is given by the following expression
$$\boldsymbol{x}(t)=\exp(\boldsymbol{A}(t-t_0))\boldsymbol{x}(t_0)+\int_{t_0}^{t}\exp(\boldsymbol{A}(t-\tau))\boldsymbol{B}f(\tau)d\tau.$$
While the integral on the right hand side is the convolution integral the expression $\exp(\boldsymbol{A}(t-t_0))$ is the matrix exponential which is defined as:
$\exp(\boldsymbol{A}(t-t_0)):=\boldsymbol{I}+\boldsymbol{A}(t-t_0)+\frac{(t-t_0)^2}{2!}\boldsymbol{A}^2+\frac{(t-t_0)^3}{3!}\boldsymbol{A}^3+\ldots .$
A: $$
my''+cy'+ky=f(t) \quad y(0)=0,\quad y'(0)=0
$$
Putting $\omega_0^2 =\frac{k}{m}$, $2\gamma=\frac{c}{m}$ and $g(t)=\frac{f(t)}{m}$ we have
$$
y''+2\gamma y'+\omega_0^2y=g(t) \quad y(0)=0,\quad y'(0)=0
$$
and applying the Laplace transform we have
$$
(s^2+2\gamma s+\omega_0^2)Y(s)=G(s)
$$
and then we can write $Y(s)=H(s)G(s)$ where
$$
H(s)=\frac{1}{s^2+2\gamma s+\omega_0^2}
$$
If we complete the square in the denominator,
$$
s^2+2\gamma s+\omega_0^2=(s+\gamma)^2 +\omega_0^2-\gamma^2
$$
we have
$$
H(s)=\frac{1}{(s+\gamma)^2 +\omega_0^2-\gamma^2}
$$
and then
$$
h(t)=\frac{\mathrm e^{-\gamma t}}{\sqrt{\omega_0^2-\gamma^2}}\sin\left(\sqrt{\omega_0^2-\gamma^2}\,t\right)
$$
and putting $\omega=\sqrt{\omega_0^2-\gamma^2}$ we can write
$$
h(t)=\frac{\mathrm e^{-\gamma t}}{\omega}\sin\left(\omega\,t\right)
$$
We now use the Convolution Property of Laplace Transforms, which states:
$$
y(t)=\mathcal L^{-1}\{H(s)G(s)\}=(h*g)(t)=\int_0^t h(t-u)g(u)\,\mathrm du
$$
that is
$$
y(t)=\frac{1}{\omega}\int_0^t \mathrm e^{-\gamma (t-u)}\sin\left(\omega\,(t-u)\right)g(u)\,\mathrm du
$$
or 
$$
y(t)=\frac{m}{\sqrt{\omega_0^2-\gamma^2}}\int_0^t \mathrm e^{-\gamma (t-u)}\sin\left(\sqrt{\omega_0^2-\gamma^2}\,(t-u)\right)f(u)\,\mathrm du
$$
