  A spectral method is developed for the direct solution of linear ordinary
differential equations with variable coefficients. The method leads to matrices
which are almost banded, and a numerical solver is presented that takes O(m^2n)
operations, where m is the number of Chebyshev points needed to resolve the
coefficients of the differential operator and n is the number of Chebyshev
coefficients needed to resolve the solution to the differential equation. We
prove stability of the method by relating it to a diagonally preconditioned
system which has a bounded condition number, in a suitable norm. For Dirichlet
boundary conditions, this implies stability in the standard 2-norm. An adaptive
QR factorization is developed to efficiently solve the resulting linear system
and automatically choose the optimal number of Chebyshev coefficients needed to
represent the solution. The resulting algorithm can efficiently and reliably
solve for solutions that require as many as a million unknowns.
