  A numerical method for solving elliptic PDEs with variable coefficients on
two-dimensional domains is presented. The method is based on high-order
composite spectral approximations and is designed for problems with smooth
solutions. The resulting system of linear equations is solved using a direct
(as opposed to iterative) solver that has optimal O(N) complexity for all
stages of the computation when applied to problems with non-oscillatory
solutions such as the Laplace and the Stokes equations. Numerical examples
demonstrate that the scheme is capable of computing solutions with relative
accuracy of $10^{-10}$ or better, even for challenging problems such as highly
oscillatory Helmholtz problems and convection-dominated convection diffusion
equations. In terms of speed, it is demonstrated that a problem with a
non-oscillatory solution that was discretized using $10^{8}$ nodes was solved
in 115 minutes on a personal work-station with two quad-core 3.3GHz CPUs. Since
the solver is direct, and the "solution operator" fits in RAM, any solves
beyond the first are very fast. In the example with $10^{8}$ unknowns, solves
require only 30 seconds.
