$N$ coupled quantum harmonic oscillators 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

(If I did not make obvious algebra errors) the elegant solution to this problem is to realize that the matrix
$$
\Lambda=\left(\begin{array}{lllr}
0&1&0\ldots&0\\
0&0&1&\ldots 0\\
\vdots&\vdots&\vdots&1\\
1&0&0\ldots&0
\end{array}\right)
$$
actually commutes with $H$ since $\Lambda$ basically maps $x_{i+1}$ to $x_i$. 
As a result, $H$ and $\Lambda$ have a common set of eigenvectors.
Since $\Lambda^n=1_{n\times n}$, the eigenvalues of $\Lambda$ satisfy $\lambda_k^n=1$ so are the $n$’th root of unity.
$$
\lambda_k=e^{2\pi i k/n}=\omega_n^k\, .
$$
The eigenvectors are then easily found to be the Fourier vectors, i.e.
$$
v_k= \frac{1}{\sqrt{n}}\left(\begin{array}{c} 
\omega_n^k\\
(\omega_n^k)^2\\
\vdots \\
1
\end{array}\right)\, .
$$
With these you can construct the matrix $U$ of eigenvectors that will diagonalise $H$.
