System of non-autonomous coupled differential equations Background: I'm trying to analytically determine the time evolution of a qubit state given an Hamiltonian $H \propto \omega \sigma_z$ and error channels for all possible erros ($\sigma_x,\sigma_y,\sigma_z$) with given error rates $\gamma_x, \gamma_y $ and $\gamma_z$. For $z$-errors I further differentiate excitation ($\propto \gamma_z n_{th}$) and de-excitation ($\propto \gamma_z (n_{th}+1)$). After costructing and solving the corresponding Lindblad master equation I end up with the following equations of motion for the states' bloch vector components (in the interaction picture):
$$
    \frac{\text{d}x(t)}{\text{d}t}  =\left[ -\gamma_z(2n_{th}+1)-4(\gamma_y + (\gamma_x-\gamma_y)\sin^2(\omega t))\right] x(t)+2\sin(2\omega t)(\gamma_x-\gamma_y)y(t)\\
    \frac{\text{d}y(t)}{\text{d}t}  =\left[ -\gamma_z(2n_{th}+1)-4(\gamma_y + (\gamma_x-\gamma_y)\cos^2(\omega t))\right] y(t)+2\sin(2\omega t)(\gamma_x-\gamma_y)x(t)\\
    \frac{\text{d}z(t)}{\text{d}t}  =-2\gamma_z((2n_{th}+1)z(t)-1)-4(\gamma_x+\gamma_y) 
$$
Solving for $z(t)$ is possible, but I haven't found a way to solve the equations for $x$ and $y$...
Problem: Solve the system of non-autonomous coupled differential equations:
$$
    \frac{\text{d}x(t)}{\text{d}t}  =\left[ b + 2a\sin^2(\omega t)\right] x(t)+a\sin(2\omega t)y(t)\\
    \frac{\text{d}y(t)}{\text{d}t}  =\left[ b + 2a\cos^2(\omega t)\right] y(t)+a\sin(2\omega t)x(t)
$$
I tried simple decoupling by substituting $u(t) = x(t) + c\cdot y(t)$ for some $c$ but this only works for the case $\omega t = \pi / 4$ ...
Generally, I have not much experience with more advanced differential equations. Is there a strategy or trick to solve them?
 A: Define $p = x+y$ and $q=x-y$. Now first add equations and then subtract them to get
$$
\frac{\text{d}x(t)}{\text{d}t} + \frac{\text{d}y(t)}{\text{d}t} =\left[ b + a\sin(2\omega t)\right]\left( x+y\right) +2a\sin^2(\omega t)(x-y)+2ay\\
\frac{\text{d}x(t)}{\text{d}t} - \frac{\text{d}y(t)}{\text{d}t}  = \left[ b + a\sin(2\omega t)\right]\left( x-y\right) +2a\sin^2(\omega t)(x+y)-2ay
$$
In other words
$$
\frac{\text{d}p(t)}{\text{d}t} =\left[ b + a\sin(2\omega t)\right]p +2a\sin^2(\omega t)q+2ay\\
\frac{\text{d}q(t)}{\text{d}t}   = \left[ b + a\sin(2\omega t)\right]q +2a\sin^2(\omega t)p-2ay
$$
Again add these two equations to get
$$
\frac{\text{d}\left(p(t) + q(t)\right)}{\text{d}t} = \left[ b + a\sin(2\omega t) + 2a\sin^2(\omega t)\right](p+q)
$$
Now introduce $\gamma = p+q$ and solve equation
$$\frac{\text{d}\gamma}{\text{d}t} = \left[ b + a\sin(2\omega t) + 2a\sin^2(\omega t)\right]\gamma$$
Like this
$$\ln \gamma = bt - \frac{a}{2\omega}\cos(2\omega t) + at - \frac{a}{4\omega}\sin(4\omega t) + c$$
where $c$ is the constant of integration. Now remember that $\gamma = p+q = (x+y) + (x-y) = 2x$ and therefore $x = \frac{(a+b)t}{2} - \frac{a}{4\omega}\cos(2\omega t) - \frac{a}{8\omega}\sin(4\omega t) + c'$. Finally replace this in one of the main equations and solve for $y(t)$.
================================ Edit =============================
You can actually decouple easily but it involves a little bit linear algebra. I was trying to show a simple way to solve it above, but now I realized you need to decouple the variables. Consider the system
$$
\begin{bmatrix} \dot{x}\\ \dot{y} \end{bmatrix} = \begin{bmatrix} b + 2a\sin^2(\omega t) & a\sin(2\omega t)\\ a\sin(2\omega t) & b + 2a\cos^2(\omega t) \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix}
$$
If you want the system to be decoupled, you need to introduce a linear map that diagonalizes the system matrix
$$
A = \begin{bmatrix} b + 2a\sin^2(\omega t) & a\sin(2\omega t)\\ a\sin(2\omega t) & b + 2a\cos^2(\omega t) \end{bmatrix}
$$
Notice that the system matrix can be written as (eigenvalue decomposition)
$$
A = \begin{bmatrix}  \cos(\omega t) &\sin(\omega t)\\ -\sin(\omega t) &\cos(\omega t) \end{bmatrix}\begin{bmatrix} b & 0\\ 0 & 2a+b\end{bmatrix}\begin{bmatrix}  \cos(\omega t) &\sin(\omega t)\\ -\sin(\omega t) &\cos(\omega t) \end{bmatrix}^T
$$
Since $\det(A) = b(b+2a)$ and $\text{tr}(A) = 2b+2a$ and you can solve the system of equations $\lambda_1 + \lambda_1 = 2b+2a$ and $\lambda_1\lambda_2 = b(b+2a)$ to find the eigen values $\lambda_1 = b$, $\lambda_2 = 2a+b$ and consequently the eigen vectors. If you want the system to always have EVD (i.e. the system matrix be diagonalizable, or equivalently, the system can be decoupled) then we must have $b(b+2a) > 0$ and $2a+2a \neq 0$, since the system matrix is symmetric already.
Therefore from EVD we find that the linear map that decouples the system is
$$
\begin{cases}x = \cos(\omega t)u + \sin(\omega t)v \\ y = -\sin(\omega t) u + \cos(\omega t)v \end{cases}
$$
Substitute $(x,y)$ in terms of $(u,v)$ in the main equations and get the decoupled system in terms of $(u,v)$. In the $u-v$ coordinate system, the equations will be decoupled.
