I am trying to solve
\begin{align*} \begin{pmatrix} \ddot{x}_{1} \\ \ddot{x}_{2} \end{pmatrix} - \begin{pmatrix} -\omega_{1}^{2} & \omega_{1}^{2} \\ \omega_{2}^{2} & -\omega_{2}^{2} \end{pmatrix}\begin{pmatrix} x_{1} \\ x_{2} \end{pmatrix} = \begin{pmatrix} a_{1} \\ a_{2} \end{pmatrix}. \end{align*}
I've been able to obtain the homogeneous solution $$ \vec{x}_{h}(t) = A\begin{pmatrix} 1 \\ 1 \end{pmatrix} + B\begin{pmatrix} t \\ t \end{pmatrix} + C\begin{pmatrix} \omega_{1}^{2}\sin\omega t \\ -\omega_{2}^{2}\sin\omega t \end{pmatrix} + D\begin{pmatrix} \omega_{1}^{2}\cos\omega t \\ -\omega_{2}^{2}\cos\omega t \end{pmatrix} $$ where $\omega = \sqrt{\omega_{1}^{2}+\omega_{2}^{2}}$, but then I am left puzzled as to how to obtain a particular solution. I tried to apply variation of parameters, but that left me with an erroneous solution (either because I made a calculation mistake or because it is inapplicable here). Does anyone have suggestions? Is there an analytic solution? If so, what is it? If not, how can we prove there isn't?
The idea is to revert to a single ODE for one of the two functions $x_1,x_2$.
Let us notice that
$$\ddot{x}_{1}-a_1 = -C( \ddot{x}_{2}-a_2) \ \text{where} \ C:=\frac{\omega_1^2}{\omega_2^2}$$
$$\overbrace{x_1+Cx_2}^{\bf{..}}=a_1+Ca_2$$
Integrating twice :
$$x_1+Cx_2=\tfrac12(a_1+Ca_2)t^2+Kt+L$$
for some constants $K,L$
It remains to plug
$$x_1=-Cx_2+\tfrac12(a_1+Ca_2)t^2+Kt+L\tag{1}$$ into the first equation
$$\ddot{x}_{1}=-\omega_1^2 x_1 +\omega_1^2 x_2 +a_1$$
to get a second order linear ordinary differential equation with constant coefficients with unknown function $x_2$ which has an exact solution. As this equation is not homogeneous (indeed, its RHS is non-zero), you have to look for a particular solution under the form of a polynomial.
Of course, once a solution for $x_2$ is found, (1) gives the expression of $x_1$.