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:** 


1. 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.
2. 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.
3. 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](https://mathoverflow.net/questions/203863/computing-gauss-legendre-quadrature-for-large-n#) 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.



Symmetric matrix formula for Gauss-Legendre quadrature
This is a particular implementation of a more general method, described in John Boyd's *Why Eigenvalues Are Roots: A Derivation of the One-Dimensional Companion Matrix for General Orthogonal Polynomials* ([restricted access](http://epubs.siam.org/doi/pdf/10.1137/1.9781611973525.appa)). Gauss-Legendre quadrature of order $n$ needs the roots of the Legendre polynomial $P\_n(x)$ and a numerical root solver must guarantee that the roots are real. By reformulating the root finding problem into an eigenvalue problem for a symmetric matrix a numerical instability for large $n$ is avoided. The specific choice of symmetric companion matrix used here is derived on page 9 of [this paper](https://cfwebprod.sandia.gov/cfdocs/CompResearch/docs/609841.pdf).


For the weights, I think the method described in the OP needs correction: It is not the $j$-th component of the eigenvector of the smallest eigenvalue that determines the weight $w\_j$, but the first component of the $j$-th eigenvector. (I looked at the code linked by the OP, and it seems that is indeed what it does.) The relation is explained on page 223 of the [Golub-Welsch paper](https://n.ethz.ch/~usoler/public/sem_work/report/literature/2004418.pdf), where all of this originates (or in more detail on page 6 of these [notes](http://www2.math.umd.edu/~mariakc/teaching/gaussian.pdf)).


