Suppose we have two quantum harmonic oscillators, with different masses $m_{1},m_{2}$ and frequencies $\omega_{1,2}$. Then we can say particles are \emph{distinguishable}, in the sense that particle $1$ has a different mass ($m_{1}$) than particle $2$ ($m_{2}$); they are not identical. Then, if the two oscillators are not coupled to each other, we can write the wavefunction for a state $|n,m\rangle$ as: \begin{eqnarray} \Psi(x_{1},x_{2}) = H_{n}(x_{1})H_{m}(x_{2}) \end{eqnarray} where $H_{n}(x)$ is the Hermite polynomial. Such wavefunction satisfies the time-independent Schrodinger equation with $$E = \hbar\omega_{1}\left(n+\frac{1}{2}\right) + \hbar\omega_{2}\left(m+\frac{1}{2}\right),$$ and represents a state where particle 1 is in the $n$th state of the oscillator, whereas particle $2$ would be in the $m$th state. I assume so far this is correct.
Now, consider the case where these two oscillators are coupled to each other, by some interacting term $V(x_{1},x_{2})$. Then, the global motion of the system is supposed to carry some entanglement, and the wavefunction does not have the simple product form as above. My question is on how could one proceed in order to find an exact representation of the wave-function in this case. Moreover, under which conditions is the system supposed to be integrable? meaning that one can compute the exact spectrum in the same way as one does for identical, indistinguishable particles (for instance in 1D, by means of Bethe ansatz methods or similars, although this is a 2-particle problem). Are there any known theorems on the integrability of such systems, for instance, considering anharmonic effects to infinite order of perturbation theory?