  In this note we provide an algorithm for computing the fractional integrals
of orthogonal polynomials, which is more stable than that using the expression
of the polynomials w.r.t. the canonical basis. This algorithm is aimed at
solving corresponding fractional differential equations. A few numerical
examples are reported.
