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*} $$