Consider the following system of coupled differential equations
$$
\dot{x}_1(t) = -x_1(t) - \cos(\omega t)x_1(t) + \cos(\omega t)x_2(t), \ x_1(0)\in\mathbb{R},\\
\dot{x}_2(t) = -\gamma x_2(t) - \cos(\omega t)x_2(t) + \cos(\omega t)x_1(t), \ x_2(0)\in\mathbb{R},
$$
where $\omega$ and $\gamma$ are positive real constants. Observe that $\bar{x}=(\bar{x}_1,\bar{x}_2)=(0,0)$ is an equilibrium of the above system.

It is almost trivial to see that if $\gamma=1$ then $\bar x$ is attractive. Indeed, in this case, we have that $x(t)=[x_1(t), x_2(t)]^\top$ can be explicitly computed as
$$
x(t) = \exp\left(\begin{bmatrix}-t &0\\ 0 & -t\end{bmatrix} + \frac{1}{\omega}\sin(\omega t)\begin{bmatrix}-1 &1\\ 1 & -1\end{bmatrix}\right)x(0),
$$
so that $x(t)\to 0$ as $t\to \infty$.

However in case $\gamma\ne 1$ proving the attractiveness of the origin is not obvious (and perhaps not even true!).

>In particular, numerical simulations seem to suggest that for $\gamma$ and $\omega$ sufficiently small (e.g. $\gamma=0.001$ and $\omega=10$) the equilibrium $\bar{x}$ is **not** attractive.

I've struggled a lot to find a way of formally proving this, with no luck. So I decided to post the problem here hoping that some of you will provide some useful suggestions or tips. Thank you!


----------


I post here the Mathematica code that I've used in my simulations:

    (* nominal values for simulation *)
    values = {gamma -> 0.001, w -> 10};
    
    equations = {
       {x1'[t], x2'[t]} == {-x1[t] - Cos[w*t]*x1[t] + Cos[w*t]*x2[t], -gamma*x2[t] - Cos[w*t]*x2[t] + Cos[w*t]*x1[t]},
       {x1[0], x2[0]} == {0.1, 0.1}};
    
    {x1t, x2t} = NDSolveValue[equations /. values, {x1[t], x2[t]}, {t, 0, 1000}];

    Plot[x1t, {t, 0, 1000}, PlotRange -> {-0.2, 0.2}]
    Plot[x2t, {t, 0, 1000}, PlotRange -> {-0.2, 0.2}]


----------

**Further remarks.** Using the change of variable approach proposed by David Hughes, i.e. $y=\begin{bmatrix}e^t & 0 \\ 0 & e^{\gamma t}\end{bmatrix}x$, we end up with the following dynamical system
$$
\dot{y}  = \cos(\omega t)\begin{bmatrix}e^t & 0 \\ 0 & e^{\gamma t}\end{bmatrix}\begin{bmatrix}-1 & 1 \\ 1 & -1\end{bmatrix} \begin{bmatrix}e^{-t} & 0 \\ 0 & e^{-\gamma t}\end{bmatrix}y = \cos(\omega t)\begin{bmatrix}-1 & e^{(1-\gamma)t} \\ e^{(\gamma-1)t} & -1\end{bmatrix}y.
$$

In view of the latter expression, we can now readily get a closed-form expression for $x(t)$, namely
\begin{align}
x(t) &= \begin{bmatrix}e^{-t} & 0 \\ 0 & e^{-\gamma t}\end{bmatrix}\exp\begin{bmatrix}-\int_0^t\cos(\omega \tau) \mathrm{d}\tau & \int_0^t\cos(\omega \tau) e^{(1-\gamma)\tau}\mathrm{d}\tau \\ \int_0^t\cos(\omega \tau)e^{(\gamma-1)\tau}\mathrm{d}\tau & -\int_0^t\cos(\omega \tau) \mathrm{d}\tau\end{bmatrix}x(0)\\
 &= \begin{bmatrix}e^{-t} & 0 \\ 0 & e^{-\gamma t}\end{bmatrix}\exp\begin{bmatrix}-\frac{1}{\omega}\sin(\omega t) & f_1(\gamma,\omega)\\ f_2(\gamma,\omega) & -\frac{1}{\omega}\sin(\omega t)\end{bmatrix}x(0),
\end{align}
where 
\begin{align}
f_1(\gamma,\omega)&:=\frac{e^{(1- \gamma)t} ((1-\gamma) \cos(\omega t) + \omega \sin(\omega t))-(1-\gamma)}{(\gamma-1)^2 + \omega^2}, \\
f_2(\gamma,\omega) &:= \frac{e^{(\gamma -1) t} ((\gamma-1) \cos(\omega t) + \omega \sin(\omega t))-(\gamma-1)}{( \gamma -1 )^2 + \omega^2}.
\end{align}

At this point, it is not clear to me how to proceed. Any ideas? 

*(A perhaps naive approach would be to compute a closed form expression of the matrix exponential above, but this seems quite nasty...)*