I want to find the wave functions of $N$ coupled quantum harmonic oscillators having the following hamiltonian:
\begin{eqnarray} H &=& \sum_{i=1}^N \left(\frac{p^2_i}{2m_i} + \frac{1}{2}m_i\omega^2 x^2_i + \frac{\kappa}{2} (x_i-x_{i+1})^2 \right)\,, \qquad x_{N+1}=0\,,\\ &=& \frac{1}{2}p^T Mp + \frac{1}{2}x^TKx\,, \end{eqnarray} where $M=\text{diag}(\frac{1}{m_1}, \cdots,\frac{1}{m_N})$ and $K$ is a real symmetric $N\times N$ matrix with positive eigenvalues, \begin{equation} K= \begin{pmatrix} k'_1& -\kappa & 0 & \cdots & 0 \\ -\kappa & k'_2& -\kappa & \ddots & \vdots \\ 0 & -\kappa & \ddots& \ddots & 0\\ \vdots & \ddots &\ddots & k'_{N-1}&-\kappa \\ 0&\cdots & 0 & -\kappa & k'_N \end{pmatrix} \end{equation} with $k'_i = m_i\omega^2+2\kappa$ but $k'_{1,N} = m_{1,N}\omega^2+\kappa$. By choosing a basis which diagonalizes the matrix $K$, the hamiltonian can be express as the sum of uncoupled harmonic oscillators hamiltonian.
As an example, consider two coupled quantum harmonic oscillators with hamiltonian \begin{equation} H = \frac{p^2_1}{2m_1} + \frac{p^2_2}{2m_2} + \frac{1}{2}m_1\omega^2 x^2_1 + \frac{1}{2}m_2 \omega^2 x^2_2 + \frac{\kappa}{2} (x_1-x_2)^2 \,. \end{equation} We make the following changes of variables (normal coordinates) \begin{eqnarray} x &=& \frac{x_1 - x_2}{\sqrt{2}} \,, \\ X &=& \frac{m_1 x_1 + m_2 x_2}{M\sqrt{2}}\,, \end{eqnarray} or equivalently, \begin{eqnarray} x_1 &=& \frac{1}{\sqrt{2}}\left(X + \frac{m_2}{M}x\right) \,, \\ x_2 &=& \frac{1}{\sqrt{2}}\left(X - \frac{m_1}{M}x\right) \,, \end{eqnarray} where $M=(m_1+m_2)/2$. Then the hamiltonian becomes \begin{equation} H = \frac{p^2_x}{2\mu} + \frac{1}{2}\mu\omega_-^2 x^2 + \frac{p^2_X}{2M} + \frac{1}{2}M\omega_+^2 X^2 \,. \end{equation} where $\displaystyle\mu = \frac{m_1m_2}{M}$ and $\omega_+^2=\omega^2$ and $\omega_-^2 = \omega^2 + 2\kappa/\mu$.
The wave functions are \begin{equation} \Psi_{mn}(x_1,x_2) = \frac{1}{\sqrt{\pi x_0X_0}}\frac{e^{-x^2/2x_0^2}}{\sqrt{m!\,2^m}}\frac{e^{-X^2/2X_0^2}}{\sqrt{n!\,2^n}}H_m\left(\frac{x}{x_0} \right)H_n\left(\frac{X}{X_0} \right) \,, \end{equation} where $x=x(x_1,x_2)$ and $X=X(x_1,x_2)$ and $\displaystyle x_0=\sqrt{\frac{\hbar}{\mu\omega_-}}$ and $\displaystyle X_0=\sqrt{\frac{\hbar}{M\omega_+}}$.
How does all this work using matrices "formalism"? And how to extend it to $N$ CQHO? Ultimately, I would like to redemonstrate (8) and (13) in http://arxiv.org/pdf/hep-th/9303048.pdf