For many applications in science and engineering, the ability to efficiently and accurately approximate solutions to elliptic PDEs dictateswhat physical phenomena can be simulated numerically. In this talk, we present a high-order accurate discretization techniquefor variable coefficient PDEs with smooth coefficients. Thetechnique comes with a nested dissection inspired direct solver that scales linearly or nearly linearly with respect to the number of unknowns. Unlike the application of nested dissection methods to classic discretization techniques, the constant prefactors do not grow with the order of the discretization. The discretization is robust even for problems with highly oscillatory solutions. For example, a problem 100 wavelengths in size can be solved to 9 digits of accuracy with 3.7 million unknowns on a desktopcomputer. The precomputation of the direct solver takes 6 minutes on a desktop computer. Then applying the computed solver takes 3 seconds. The recent application of the algorithm to inverse media scattering also will be presented.