  A numerical method for variable coefficient elliptic problems on two
dimensional domains is described. The method is based on high-order spectral
approximations and is designed for problems with smooth solutions. The
resulting system of linear equations is solved using a direct solver with
$O(N^{1.5})$ complexity for the pre-computation and $O(N \log N)$ complexity
for the solve. The fact that the solver is direct is a principal feature of the
scheme, and makes it particularly well suited to solving problems for which
iterative solvers struggle; in particular for problems with highly oscillatory
solutions. This note is intended as a tutorial description of the scheme, and
draws heavily on previously published material.
