We consider the problem of numerically generating the recursion coefficients of orthogonal polynomials, given an arbitrary weight distribution of either discrete, continuous, or mixed type. We discuss two classical methods, respectively due to Stieltjes and Chebyshev, and modern implementations of them, placing particular emphasis on their numerical stability properties. The latter are being studied by analyzing the numerical condition of appropriate finite-dimensional maps. A number of examples are given to illustrate the strengths and weaknesses of the various methods and to test the theory developed for them.