  A fast algorithm (linear in the degrees of freedom) for the solution of
linear variable-coefficient rational-order fractional integral and differential
equations is described. The approach is related to the ultraspherical method
for ordinary differential equations [Olver & Townsend 2013], and involves
constructing two different bases, one for the domain of the operator and one
for the range of the operator. The bases are constructed from direct sums of
suitably weighted ultraspherical or Jacobi polynomial expansions, for which
explicit representations of fractional integrals and derivatives are known, and
are carefully chosen so that the resulting operators are banded or
almost-banded. Geometric convergence is demonstrated for numerous model
problems when the variable coefficients and right-hand side are sufficiently
smooth.
