Gauss-Legendre Quadrature, computation of the abscissas and weights I would like to write a program to calculate abscissas  and weights of Gauss-Legendre Quadrature. 
I found the following source
http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature.
I don't understand why the first guess $x_0$ for the $i$-th root of a $n$-order polynomial $P_n$  is defined as
$x_0=\cos\left(\pi\frac{i-\frac{1}{4}}{n+\frac{1}{2}}\right)$
Could someone explain this to me?
 A: H. Bruns proved in 1881 that the arccosine of the $i$th root of $P_n$ lies in the interval 
$$
\left[ \frac{2i-1}{2n+1}\pi, \ \frac{2i}{2n+1}\pi\right]
$$
(When working with Legendre polynomials, it is often convenient to compose them with cosine, i.e., look at $P_n(\cos\theta)$ instead of $P_n(x)$.) 
So it's natural to take the cosine of the midpoint of this interval as the first guess for the root. One source for the above is Szegő's paper mentioned by uranix; another one is Szegő's book Orthogonal Polynomials, Theorem 6.21.2.
That said, I think that using a root search procedure for each root individually is not particularly efficient, as this does not take advantage of the special structure of the polynomials $P_n$. The recurrence property of $P_n$ leads to a family of tridiagonal matrices for which $P_n$ is (a multiple of) the characteristic polynomial. It's easier to find the eigenvalues of a tridiagonal matrix than it is to find the roots of a generic polynomial: see the paper Calculation of Gauss quadrature rules by G. H. Golub and J. H. Welsch. 
Spoiler: the paper Is Gauss Quadrature Better
than Clenshaw–Curtis? by Lloyd N. Trefethen has the following 7-line implementation of Gaussian quadrature on $[-1,1]$ in Matlab. 
function I = gauss(f,n)               % (n+1)-pt Gauss quadrature of f
beta = .5./sqrt(1-(2*(1:n)).ˆ(-2));   % 3-term recurrence coeffs
T = diag(beta,1) + diag(beta,-1);     % Jacobi matrix
[V,D] = eig(T);                       % eigenvalue decomposition
x = diag(D); [x,i] = sort(x);         % nodes (= Legendre points)
w = 2*V(1,i).ˆ2;                      % weights
I = w*feval(f,x);                     % the integral

