Let $H(r, R)$ be a hamiltonian of a system in center-of-mass and relative coordinates, being $r = r_1 - r_2$ and $$R = \frac{m_1r_1 + m_2 r_2}{m_1 + m_2}$$ Consider the Hamiltonian to contain coupling terms between $r$ and $R$ as expressed in the next equation $$H(r,R) = H_{rm}(r) + H_{CM}(R) + rR.$$ The question is how can we find the complete basis of wavefunctions that diagonalize the total Hamiltonian?
My suggestion is to first diagonalize the relative and CM components, so we can build a new basis of the form $$\Psi_{n,m}(r,R) = \psi_n(r) \phi_m(R)$$ Therefore we can diagonalize the coupling term as $$\left< \Psi_{n,m}(r,R) \right| rR \left|\Psi_{n,m}(r,R)\right>$$ Is it right or are there any other way to find the complete solution?