A prototypical model for condensed-phase chemical reactions is the spin-boson model: a two-level spin coupled to a continuum of vibrations (bosons). The spin-boson Hamiltonian is
\[ H = \sigma_x + \sigma_z \otimes \int_0^{x_\mathrm{max}} h(x) \rho(x)(a_{x}^\dagger + a_{x})\mathrm{d}x + \int_0^{x_\mathrm{max}} x\rho(x) a_{x}^\dagger a_{x} \mathrm{d}x \]
Assuming $\rho(x)=1$, to discretize the continuum of bosons, we approximate the integral $\int_0^{x_\mathrm{max}} h(x) (a_{x}^\dagger + a_{x})\mathrm{d}x$ by
\[ \sum_{n=1}^N \left\{ h(x_n) (a_{x_n}^\dagger + a_{x_n})W_n \right\} \]
where
\[ W_n = \int_0^{x_\mathrm{max}} \frac{\prod_{m\neq n}(x-x_m)}{\prod_{m\neq n}(x_n-x_m)} \mathrm{d}x \]
The integrand of $W_n$ is a Lagrange interpolating polynomial. $W_n$ is the weight for the Gaussian quadrature. It has been shown1 2 that the approximation (when N<\infty) is the best if the nodes $\{x_m\}$ are the roots of the N-th order Legendre polynomials over the integration interval $[0,x_\mathrm{max}]$.
Similarly, the integral $\int_0^{x_\mathrm{max}} x\rho(x) a_{x}^\dagger a_{x} \mathrm{d}x$ can be approximated by
\[ \sum_{n=1}^N \left\{ x_n (a_{x_n}^\dagger a_{x_n})W_n \right\}. \]
To obtain the nodes, which are the roots of the N-th order Legendre polynomials, and the weights associated with the roots, it is helpful to use the recurrence relation of the Legendre polynomials. The recurrence relation for the monic Legendre polynomials over the interval $[-1,1]$ is
\[ p_{n+1}(x) = x p_{n}(x) - \frac{n^2}{4n^2 -1} p_{n-1}(x) \]
Replacing the variable $x$ with $x=\frac{2(y-a)}{b-a} -1$, and defining $\tilde{p}_n(y)=p_n\left(\frac{2(y-a)}{b-a} -1\right)$, the following recurrence is valid over the interval $[a,b]$.
\[ \tilde{p}_{n+1}(y) = \left[\frac{2(y-a)}{b-a} -1 \right] \tilde{p}_{n}(y) - \frac{n^2}{4n^2 -1} \tilde{p}_{n-1}(y) \]
The polynomial $\tilde{p}_{n}(y)$ is not monic. To make it monic, we divide it by $(b-a)^n/2^n$ and obtain a new monic polynomial on the interval $[a,b]$:
\[ p_{n}(y) = \left[ y- \frac{a+b}{2}\right]p_n(y) - \left(\frac{b-a}{2}\right)^2 \frac{n^2}{4n^2 -1} p_{n-1}(y) \]
$p_n(y)$ is, in fact, the product $\prod_{k=1}^n(y-y_k)$ where $y_k$ is the root of $p_n(y)$.
Diagonalizing the tridiagonal matrix formed by the recurrence coefficients $a_n=\frac{a+b}{2}$ and $b_n=\left(\frac{b-a}{2}\right)^2 \frac{n^2}{4n^2 -1}$ produces the nodes and the weights. The matrix is
\[ \left[\begin{array}{ccccc} a_{0} & \sqrt{b_{1}} & & & \\ \sqrt{b_{1}} & a_{1} & \sqrt{b_{2}} & & \\ & \vdots & \vdots & & \\ & & \sqrt{b_{N-2}} & a_{N-2} & \sqrt{b_{N-1}} \\ & & & \sqrt{b_{N-1}} & a_{N-1} \end{array}\right] \]
The roots are the eigenvalues of this matrix. Supposing the eigenvectors are $\mathbf{v}_n$, the weights $W_n$ are $(b-a)v_{n,1}^2$.