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). 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).