Implement Golub-Welsch Algorithm for Numerical Quadrature Rule Generation

[wavefunction91]

For some weight functions, (e.g. $\ln^2$ of Gill and Chien), there are not simple methods for the computation of quadrature nodes and weights. The Golub-Welsch algorithm provides a systematic way of numerically calculating these quantities for arbitrary weight functions (given the ability to compute moments of specified order).

It requires the solution of a Hankel-type tridiagonal eigenvalue problem, which tends to be numerically unstable for large orders, but having it as an option would lower the barrier to entry for testing out new quadratures and an important validation tool moving forward.

Requires:

• Integration of LAPACK++ to solve the EVP
• Tabulate coefficients for the classical polynomials