I am trying to solve the following system of differential equations
$$\begin{bmatrix}\frac{dP_-}{dt}\\\frac{dP_+}{dt}\end{bmatrix}=\begin{bmatrix}-Γ_+&Γ_-\\Γ_+&-Γ_-\end{bmatrix}\begin{bmatrix}P_-\\P_+\end{bmatrix},$$
where $Γ_+=Γ^cn_p$, $Γ_-=Γ^c(1+n_p)$ and
$$n_p = \left[\mathrm{exp \left( \frac{ \sqrt{\Delta_c^2+(\epsilon_0+\delta \epsilon \sin(\omega t))^2}}{k_B T}\right)}-1\right]^{-1}$$
holds and $P_-+P_+ = 1$. You can see that with these $n_p$ factors the coefficients become time-dependent. I need to find the steady-state solution for $P_-$ and $P_+$ to first order in $\delta \epsilon / \Delta_c $. For more details and physical motiviation see this paper.
I already realized that to find the solution I only need one of the equations, but I don't know how should I find the actual result mentioned in the paper.
Currently I am trying with the following idea: use the first equation of the system and the constraint that $P_- + P_+ = 1$, to write
\begin{equation}
\frac{dP_-}{dt} = -(\Gamma_+ + \Gamma_-)P_- + \Gamma_-
\end{equation}
and solve this. I think this is a called a non-homogeneous differential equation, which means that I should find the general solution of the homogeneous part and a particular solution for the non-homogeneous part. Actually the solution should be independent of the initial condition which gives me a bit more confusion.
Maybe I should post this question on the physics StackExchange, sorry for the
inconvenience.
Best Answer
So after some calculation I found that there is a well-defined expression for solving differential equations like this. This expression can be found by the method of constant variation. However, in my case this expression cannot be evaluated exactly (the integrals cannot be done), this is why we need to Taylor expand the integrands to first order in the small $\delta \epsilon / \Delta_c$ quantity. Then the solution can be found by evaluating the integrals.