To understand how the weights and the mesh points are generated, we define first a polynomial of degree 2N-1 (since we have 2N variables at hand, the mesh points and weights for N points). This polynomial can be represented through polynomial division by \begin{equation*} P_{2N-1}(x)=L_N(x)P_{N-1}(x)+Q_{N-1}(x), \end{equation*} where P_{N-1}(x) and Q_{N-1}(x) are some polynomials of degree N-1 or less. The function L_N(x) is a Legendre polynomial of order N .
Recall that we wanted to approximate an arbitrary function f(x) with a polynomial P_{2N-1} in order to evaluate \begin{equation*} \int_{-1}^1f(x)dx\approx \int_{-1}^1P_{2N-1}(x)dx. \end{equation*}