One analytical model of the marriage relationship — Romeo & Juliet — includes the constant term before marriage (uninfluenced state), call it, $x_0$ for one spouse and $y_0$ for the other spouse. Let $r_1$ and $r_2$ be the exponential rates of 'decay' for love or hate. One spouse may tend toward the uninfluenced state $x_0$ so the ODE is $\dot x = r_1 (x_0 - x)$ and $\dot y = r_1 (y_0 - y)$. Now add in the influence function $I(x)$ and $I(y)$ where $I(x) = a_1 y$ and $I(y) = a_2 x$. Therefore the coupled ODE is
$$ \begin{aligned} \dot x &= r_1 (x_0 - x) + a_1 y \\ \dot y &= r_2 (y_0 - y) + a_2 x \end{aligned} $$
How to solve this system of coupled linear ODEs? 1) Solve for $y$ in Eq (1) then sub into Eq. (2)??
First, that system of equations IS linear.
I would start by differentiating each equation again: $x''= -r_1x'+ a_1y'$ $y''= -r_1y'+ a_2x'$
Replace the $y'$ in the first equation by $y'= r_2(y_0- y)+ a_2x$: $x''= -r_1x'+ a_1(r_2(y_0- y)+ a_2x)= (a_1a_2- r_1)x- a_1r_2y_0+ a_1r_2y$
Finally replace that remaining y with $(x'+ r_1x)/a_1$: $x''= (a_1a_2- r_1)x- a_1r_2y_0+ r_2x'+ r_1r_2x$
So that we have the single linear differential equation in x with constant coefficients:
$x''- r_2x'- (r_1r_2+ a_1a_2- r_1)x= a_1r_2y_0$.