Solving a set of coupled first order differential equations I have the following set of coupled equations (which results from consideration of the motion of a particle in a one dimensional non-autonomous harmonic potential):
$$\dot{q}(t)=\frac{p(t)}{m}e^{-\gamma t}\\ ~~~~~~~\dot{p}(t)=-m\omega^2q(t)e^{\gamma t} $$
I am asked to give a general solution of the equation of motion in the case $0<\gamma<\omega$. A dot represents a total time derivative.
I thought I might be able to rewrite it like this:
$$\frac{d}{dt}\begin{bmatrix}p\\q\end{bmatrix}=\begin{bmatrix}0&&\frac{1}{m}e^{-\gamma t}\\-m\omega^2e^{\gamma t}&&0\end{bmatrix}\begin{bmatrix}p\\q\end{bmatrix}\equiv A\begin{bmatrix}p\\q\end{bmatrix}$$
Which has solutions $$\begin{bmatrix}p(t)\\q(t)\end{bmatrix}=e^{tA}\begin{bmatrix}q_0\\p_0\end{bmatrix}$$Where $\{q_0,p_0\}$ are determined by the initial conditions. However I am not sure if this is true for a time dependant matrix A. I suspect that this is not true so I abandoned this method.  
Next I tried guessing a solution of the form $q(t)=e^{\lambda t}$ and plugged this into the two differential equations which told me the two possible values of $\lambda$. Then the general solution would be linear combinations of those two guessed solutions, i.e:
$$q(t)=C_1e^{\lambda_+t}+C_2e^{\lambda_-t}\\=C_1\exp{[\frac{1}{2}t(-\gamma +\sqrt{\gamma^2-4\omega^2})]}+C_2\exp{[\frac{1}{2}t(-\gamma -\sqrt{\gamma^2-4\omega^2})]}$$
But Since we are considering the case $0<\gamma<\omega$, the argument of the radical is negative and so the $\lambda_{\pm}$ are complex. So I could rewrite this as:
$$q(t)=e^{-\frac{1}{2}\gamma t}(C_1e^{i\alpha t}+C_2e^{-i\alpha t})$$
Where $\alpha=\frac{1}{2}\sqrt{4\omega^2-\gamma^2}>0$. But I wouldn't know how to physically interpret this solution (since it is complex and not real) in all cases except $(C_1,C_2)=(1,-1)$ which would allow us to rewrite it as $\sin{\alpha t}$ using Euler's formula. Which would mean we are only allowed (very) specific initial conditions. So either my method is incorrect or I am interpreting it incorrectly. Which one is it? If it is my method could you point me in the right direction? (Also I would like to know if my first method is valid, i.e can you use the matrix exponential in the non-autonomous case). 
And if anyone was interested, the original Lagrangian of the system was
$$\mathcal{L}=\frac{1}{2}me^{\gamma t}(\dot{q}^2-\omega^2q^2),$$
For a constant $\omega$ and $\gamma$.
 A: Well, we can use Laplace transform:
$$\mathscr{L}_t\left[-\text{m}\cdot\omega^2\cdot\exp\left(\gamma\cdot t\right)\cdot\text{q}\left(t\right)\right]_{\left(\text{s}\right)}=-\text{m}\cdot\omega^2\cdot\mathscr{L}_t\left[\exp\left(\gamma\cdot t\right)\cdot\text{q}\left(t\right)\right]_{\left(\text{s}\right)}=$$
$$-\text{m}\cdot\omega^2\cdot\mathscr{L}_t\left[\text{q}\left(t\right)\right]_{\left(\text{s}-\gamma\right)}=-\text{m}\cdot\omega^2\cdot\text{Q}\left(\text{s}-\gamma\right)\tag1$$
And:
$$\mathscr{L}_t\left[\text{p}'\left(t\right)\right]_{\left(\text{s}\right)}=\text{s}\cdot\text{P}\left(\text{s}\right)-\text{p}\left(0\right)\tag2$$
Now, for the first equation we can use the same:
$$\mathscr{L}_t\left[\frac{\text{p}\left(t\right)}{\text{m}}\cdot\exp\left(-\gamma\cdot t\right)\right]_{\left(\text{s}\right)}=\frac{1}{\text{m}}\cdot\mathscr{L}_t\left[\text{p}\left(t\right)\right]_{\left(\text{s}+\gamma\right)}=\frac{1}{\text{m}}\cdot\text{P}\left(\text{s}+\gamma\right)\tag3$$
And:
$$\mathscr{L}_t\left[\text{q}'\left(t\right)\right]_{\left(\text{s}\right)}=\text{s}\cdot\text{Q}\left(\text{s}\right)-\text{q}\left(0\right)\tag4$$
So, we get:
$$
\begin{cases}
\text{s}\cdot\text{Q}\left(\text{s}\right)-\text{q}\left(0\right)=\frac{1}{\text{m}}\cdot\text{P}\left(\text{s}+\gamma\right)\\
\\
\text{s}\cdot\text{P}\left(\text{s}\right)-\text{p}\left(0\right)=-\text{m}\cdot\omega^2\cdot\text{Q}\left(\text{s}-\gamma\right)
\end{cases}\tag5
$$
So, we get:
$$\text{s}\cdot\frac{\text{p}\left(0\right)-\text{s}\cdot\text{P}\left(\text{s}+\gamma\right)}{\text{m}\cdot\omega^2}-\text{q}\left(0\right)=\frac{1}{\text{m}}\cdot\text{P}\left(\text{s}+\gamma\right)\space\Longleftrightarrow\space$$
$$\text{P}\left(\text{s}+\gamma\right)=\frac{\text{s}\cdot\text{p}\left(0\right)-\text{m}\cdot\omega^2\cdot\text{q}\left(0\right)}{\text{s}^2+\omega^2}\tag6$$
A: $$\begin{cases}\frac{dq}{dt}=\frac{p(t)}{m}e^{-\gamma t} \quad\to\quad p(t)=me^{\gamma t}\frac{dq}{dt}\\ \frac{dp}{dt}=-m\omega^2q(t)e^{\gamma t}\end{cases} $$
$$\frac{dp}{dt}=me^{\gamma t}\left(\gamma\frac{dq}{dt}+\frac{d^2q}{dt^2}\right)=-m\omega^2q(t)e^{\gamma t} $$
$$\frac{d^2q}{dt^2}+\gamma\frac{dq}{dt}+\omega^2q=0 $$
This is a second order linear ODE easy to solve.
