A procedure is proposed in order to expand $w={\prod }_{j=1}^N{P}_{i_j}\left(x\right)={\sum }_{k=0}^M{L}_k{P}_k\left(x\right)$ where $P_i\left(x\right)$ belongs to aclassical orthogonal polynomial sequence (Jacobi, Bessel, Laguerre and Hermite) ($M={\sum }_{j=1}^N{i}_j$). We first derive a linear differential equation of order $2^N$ satisfied by w, fromwhich we deduce a recurrence relation in k for the linearizationcoefficients $L_k$. We develop in detail the two cases ${\left[P_i\left(x\right)\right]}^N$, $P_i\left(x\right)P_j\left(x\right)P_k\left(x\right)$ and give the recurrencerelation in some cases (N=3,4), when the polynomials $P_i\left(x\right)$are monic Hermite orthogonal polynomials.

Let $P_k$ be any sequence of classical orthogonal polynomials of a discrete variable. We give explicitly a recurrence relation (in k) for the coefficients in $P_i{P}_j={\sum }_{k}c(i,j,k)P_k$, in terms of the coefficients σ and τ of the Pearson equation satisfied by the weight function ϱ, and the coefficients of the three-term recurrence relation and of two structure relations obeyed by $P_k$.