  We present a spectral method for one-sided linear fractional integral
equations on a closed interval that achieves exponentially fast convergence for
a variety of equations, including ones with irrational order, multiple
fractional orders, non-trivial variable coefficients, and initial-boundary
conditions. The method uses an orthogonal basis that we refer to as Jacobi
fractional polynomials, which are obtained from an appropriate change of
variable in weighted classical Jacobi polynomials. New algorithms for building
the matrices used to represent fractional integration operators are presented
and compared. Even though these algorithms are unstable and require the use of
high-precision computations, the spectral method nonetheless yields
well-conditioned linear systems and is therefore stable and efficient. For
time-fractional heat and wave equations, we show that our method (which is not
sparse but uses an orthogonal basis) outperforms a sparse spectral method
(which uses a basis that is not orthogonal) due to its superior stability.
