The classical orthogonal polynomials are given as the polynomial solutions p n ( x) of the differential equation σ(x)y″(x)+τ(x)y ′(x)+λ ny(x)=0, where σ( x) is a polynomial of at most second degree and τ( x) is a polynomial of first degree. In this paper a general method to express the coefficients A n, B n and C n of the recurrence equation p n+1(x)=(A nx+B n)p n(x)−C np n−1(x) in terms of the given polynomials σ( x) and τ( x) is used to present an algorithm to determine the classical orthogonal polynomial solutions of any given holonomic three-term recurrence equation, i.e., a homogeneous linear three-term recurrence equation with polynomial coefficients. In a similar way, classical discrete orthogonal polynomial solutions of holonomic three-term recurrence equations can be determined by considering their corresponding difference equation σ(x)Δ∇y(x)+τ(x)Δy(x)+λ ny(x)=0, where Δy( x)= y( x+1)− y( x) and ∇ y( x)= y( x)− y( x−1) denote the forward and backward difference operators, respectively, and a similar approach applies to classical q-orthogonal polynomials, being solutions of the q-difference equation σ(x)D qD 1/qy(x)+τ(x)D qy(x)+λ q,ny(x)=0, where D qf(x)= f(qx)−f(x) (q−1)x , q≠1, denotes the q-difference operator.