Solution (source code)

= Solution

Introduce the fixed endpoint values $x_0=x_{N+1}=0$. The <Lagrangian> is
$$
\boxed{L=\frac m2\sum_{j=1}^N\dot x_j^2-\frac k2\sum_{j=0}^N(x_{j+1}-x_j)^2.}
$$
Thus $T=mI_N$, $V_{ij}=k(2\delta_{ij}-\delta_{i,j+1}-\delta_{i,j-1})$, and
$$
\boxed{V-\omega^2T=\begin{pmatrix}d&-k&0&\cdots&0\\-k&d&-k&\ddots&\vdots\\0&-k&d&\ddots&0\\\vdots&\ddots&\ddots&\ddots&-k\\0&\cdots&0&-k&d\end{pmatrix},\qquad d=2k-m\omega^2.}
$$
This is a symmetric tridiagonal matrix, reflecting that only nearest-neighbour masses are coupled.