Numerical Computation of Orthogonal Polynomials Recurrence Relations Background and notations: Given an interval $I\subseteq \mathbb{R}$ and a continuous finite measure $d\mu = w(x)dx$, and denote $p_n(x)$ the orthogonal polynomials with respect to $d\mu$. We have the following recurrence relation $$ p_j (x) = (a_j x +b_j)p_{j-1}(x) +c_jp_{j-2}(x) \, ,\quad \forall j\geq 1\, \quad p_0(x) \equiv 1,\,p_{-1}(x)\equiv 0 \, ,$$ where the constants are determined by the measure.
In the famous Golub-Welsch paper, section 4, they give a numerical method to calculate the constants. However, it requires the numerical value of $\int\limits_{I}x^\ell w(x)\,dx$ for all non-negative integers $\ell$.
The problem is that for numerical integration we usually need a quadrature formula, for which weights we need the recurrence relation (see the other sections of the same paper, for example).
Question 1: Is there a way to compute $\int\limits_{I}x^\ell w(x)\,dx$ without any quadrature formula?
Question 2: Is there a way to compute the recurrence constants without evaluating integrals?
Remark: We can always evaluate $\int\limits_{I}x^\ell w(x)\,dx$ using the Gauss Legendre quadrature. This means that we need only the well-known recurrence for the Lebesgue measure to compute these integrals for all other continuous measures. I'm looking for something else, though.
 A: Computation of the coefficients of the recurrence relations of orthogonal polynomials is studied in details in the standart reference:
Gautschi, W., Orthogonal polynomials: computation and approximation. 
Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2004
In particular, discretization methods are described in Section 2.2.
They consist in the following :
1) approach the given measure $d\mu$ by a discrete $N$-point measure $d\mu_{N}$, 
2) compute the recurrence coefficients of $d\mu_{N}$ and let $N$ go to infinity.
The three main issues, namely, appropriate choice of discretization, computation of recurrence coefficients of discrete measures, and convergence as $N$ tends to infinity, are discussed in Sections 2.2.1--2.2.4. Nontrivial examples are also given.
The following recent book should also be of interest :
Gautschi, W., Orthogonal polynomials in MATLAB. Exercises and solutions. Software, Environments, and Tools, 26. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016.
A: I'm not sure if this is what you are looking for, but here we go. From orthogonality, you get 
$$\int_Ip_n(x)p_m(x)\,d\mu(x)=h_n\delta_{m,n}.$$
Suppose now that the leading coefficient of $p_n(x)$ is $k_n$. Then,
$$a_n=\frac{k_{n+1}}{k_n} \qquad \text{and} \qquad -c_{n+1}=\frac{a_{n+1}h_{n+1}}{a_nh_n}.$$
