Numerical integration is among the most common tasks in astrophysics. Simple formulae, like trapezoidal, Newton-Cotes or Simpson, are often not enough accurate. Gauss quadrature provides a more accurate solution, but its implementation is a bit more difficult. Here I will demonstrate the implementation of that algorithm in IDL. The algorithm includes finding the coefficients of Legendre polynomials and their zeros. Both tasks are explicitly solved in the following as well.
But let's go step by step. Our task is to solve numerically the integral:
$$I = \int_a^b f(x)\, \mathrm{d} x.$$
