Symmetric matrix formula for Gauss-Legendre quadrature 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.
 A: 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). 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.
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, where all of this originates (or in more detail on page 6 of these notes).
A: You are describing the Golub-Welsch algorithm. As described, it costs O(n^2) operations to compute n-point Gauss-Legendre quadrature rule.
http://www.ams.org/journals/mcom/1969-23-106/S0025-5718-69-99647-1/
If you are interested in large Gauss-Legendre quadrature rule, then better (faster and more accurate) algorithms now exist.  See your previous question: 
Computing Gauss Legendre Quadrature for Large N
