I am having trouble finding the eigenvalues for the Hamiltonian
$$ H = \frac{P_1^2}{2M} + \frac{P_2^2}{2m} + \frac{K}{2}x_1^2 + \frac{k}{2}(x_1 - x_2)^2$$
Even though I can find a basis where the $x_1$ and $x_2$ coordinates are decoupled I then get products of momenta in my new Hamiltonian. I think this problem should be transformable into two decoupled oscillators. Am I wrong about that?
Structurally, it looks like this:
Then you're done - you've got two decoupled oscillators whose position and momentum variables are canonically conjugate to each other and given as explicit linear combinations of the old ones.
The Hamiltonian is a quadratic form both in the momenta and in the coordinates. You can decouple the problem by using the principal axis theorem on the coordinates to obtain generalized coordinates. For this, first a scaling of the coordinates to equalize the coefficients of the kinetic energy is necessary. The whole approach corresponds to "diagonalizing" both the matrix of the potential energy and of the kinetic energy by finding a suitable linear transformation consisting of a scaling and an subsequent orthogonal transformation. The latter involves finding the eigenvalues and eigenvectors of the potential energy matrix.
Note: Following the comment by @Emilio Pisanty, I edited the previous cursory answer and removed some ambiguities pointed out by him. This is directed at physicists who have absolved only a standard curriculum, and thus probably have not indulged in "symplectomorphisms".
(1) The problem corresponds to the typical small oscillation problem of coupled oscillators which is treated in most classical mechanics courses, usually in the framework of the Lagrange function, which can easily be translated to the corresponding Hamilton function. The decoupling of the given Hamilton function into a sum of uncoupled linear oscillator Hamilton functions is the problem of finding the "normal coordinates" for these oscillators.
This is described in detail, e.g. in H. Goldstein's book, Classical Mechanics, Chapter 6. Oscillations, especially, section 6.2 The Eigenvalue Equation and the Principal Axis Transformation. In this Harvard lecture, you'll find a shorter exposition of Goldstein's approach. It is completely sufficient to perform the transformations on the classical Hamiltonian (or Lagrange function). Once you have the decoupled Hamiltonian, you can quantize it. Goldstein designates the resulting necessary total transformation as "similarity transformation" to realize the "principal axis transformation".
(2) In order to decouple the oscillators, you have to simultaneously diagonalize the symmetric quadratic form (matrix) of the potential energy and the symmetric quadratic form (matrix) of the kinetic energy of the Lagrangian $L=T-V$ (expressed in velocities) by a suitable linear coordinate transformation. If the masses in the kinetic energy term were equal, you could do this with a single orthogonal coordinate transformation, usually called principal axis transformation.
(3) But here the masses are not equal. Therefore, you have first to scale the coordinates so that the coefficients of $v_1^2$ and $v_2^2$ become equal. This means, that the kinetic energy matrix is not only diagonal but has also the same coefficients. Then, after a subsequent orthogonal transformation, the quadratic form of the kinetic energy will stay a sum of the squares of the velocity with the same factor. This is described in an easily understandable way (directly applicable to the present 2-D case) in P. K. Aravind, Geometrical interpretation of the simultaneous diagonalization of two quadratic forms, American Journal of Physics 57, 309 (1989). I don't know whether there is any no-pay option to get this beautiful paper.
(4) Now you perform the orthogonal transformation ( this case, rotation in the $x_1$, $x_2$ plane) to diagonalize the symmetric matrix of the potential energy. The matrix of the kinetic energy stays diagonalized. Thus you end up with both matrices diagonalized so that the Lagrangian becomes a sum of separate oscillator Lagrangians.
(5) You can find the resulting Hamiltonian $H$ in the new generalized coordinates from the Lagrangian $L$ in the new generalized coordinates by the known relation $$H= \sum_i p_i \dot x_i-L$$ and $$p_i= \frac {\partial L}{\partial \dot x_i}$$ This ensures that your new generalized momenta $p_i$ are the canonical conjugates of the new generalized coordinates $x_i$.
(6) If you need the QM (Schrödinger) Hamiltonian, you quantize the classical one by substituting $p \to -i\hbar \nabla$.
Another method to solve the problems of this kind is to consider operators of the form $$ \hat{a} = c x_1 + d x_2 + e P_1 + f P_2. $$ Here $c$, $d$, $e$ and $f$ are numerical coefficients. If you find solutions of the equation $$ \left[ \hat{a}, H \right] = \varepsilon \hat{a}, $$ you will get creation and annihilation operators and energies of quanta.
More generally, given a semipositive definite quadratic real Hamiltonian $$H=\frac{1}{2}z^I H_{IJ} z^J~\geq~ 0 \tag{1}$$ in $2n$ canonical coordinates $$(z^1, \ldots z^{2n})~=~(q^1,\ldots, q^n,p_1, \ldots, p_n),\tag{2}$$ one may show that there always exists a real & linear symplectic transformation $$Z~=~ S z, \qquad S\in Sp(2n,\mathbb{R}),\tag{3}$$ that brings the Hamiltonian on diagonal form, cf. Ref. 1.
There are related Phys.SE posts here & here.
References: