GAUSS-LEGENDRE QUADRATURE.

CASE: n = 2

The roots of the Legendre polynomial P_n are the quadrature nodes.

r := Roots[LegendreP[2,x]==0,x]

r = { r[[1]][[2]],r[[2]][[2]]}

Construct kernel matrix using all Legendre polynomials up to P_{n-1} at the quadrature nodes.

K={ {LegendreP[0,r[[1]]], LegendreP[0,r[[2]]]}, {LegendreP[1,r[[1]]], LegendreP[1,r[[2]]]}};

MatrixForm[K]

Scaling factor for sampling matrix which corrects for the scaling factor of the first basis function coefficient, which, up to this constant, is the sought integral.

n :=Integrate[LegendreP[0,x],{x,-1,1}]

Sampling matrix, with scaling factor. The first column are the weights of Gauss-Legendre quadrature (compare for example with http://en.wikipedia.org/wiki/Gaussian_quadrature).

S := n* Inverse[K]

MatrixForm[S]

CASE: n = 3

r := Roots[LegendreP[3,x]==0,x]

r = { r[[1]][[2]],r[[2]][[2]], r[[3]][[2]]}

K={ {LegendreP[0,r[[1]]], LegendreP[0,r[[2]]],LegendreP[0,r[[3]]]},
{LegendreP[1,r[[1]]], LegendreP[1,r[[2]]],LegendreP[1,r[[3]]]},
{LegendreP[2,r[[1]]], LegendreP[2,r[[2]]],LegendreP[2,r[[3]]]}};

MatrixForm[K]

n := Integrate[LegendreP[0,x],{x,-1,1}];

S:= n * Inverse[K];

MatrixForm[S]