  This paper presents a numerical method for variable coefficient elliptic PDEs
with mostly smooth solutions on two dimensional domains. The PDE is discretized
via a multi-domain spectral collocation method of high local order (order 30
and higher have been tested and work well). Local mesh refinement results in
highly accurate solutions even in the presence of local irregular behavior due
to corner singularities, localized loads, etc. The system of linear equations
attained upon discretization is solved using a direct (as opposed to iterative)
solver with $O(N^{1.5})$ complexity for the factorization stage and $O(N \log
N)$ complexity for the solve. The scheme is ideally suited for executing the
elliptic solve required when parabolic problems are discretized via
time-implicit techniques. In situations where the geometry remains unchanged
between time-steps, very fast execution speeds are obtained since the solution
operator for each implicit solve can be pre-computed.
