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 [27], 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.