While searching the web, I came across the following algorithm for the Gauss-Legendre quadrature. I wasn't able to find a reference or a proof of my own as to why it works. I'll present it, and the question is - Why does it work?

**The Method:**

We construct the

**Matrix**$A \in M_{N\times N} (\mathbb{R} )$ by $\forall 1\leq j\leq N-1 \, A_{j,j+1} = A_{j+1,j} = \frac{j}{\sqrt{4j^2 -1}}$, and all other entries are zero.For some reason, the Legendre polynomial of order $N$ is the characteristic polynomial $p_N(x) = {\rm det}(x{\rm Id - A})$. Therefore, the eigenvalues of $A$ are the

**abscissas**of the quadrature.For the

**weights**, we take $v^j$ to be the eigenvector of the $j-th$ smallest (negative) eigenvalue, and we have that $w_j =2(v^j _1) ^2, \, \forall 1\leq j \leq N$, with $v^j_1$ denoting the first coordinate.

I don't understand this method, nor can I find a proper reference to it in standard Numerical Integration books (Davis and Rebinowitz) or Orthogonal Polynomials books (Szego).

*Meta - Note:* in my old question I asked for *a* stable and efficient algorithm for the Gauß-Legendre quadrature, but none of the proposed methods was the one I've found. This is why I've opened a new question.