  A numerical scheme is presented for approximating fractional order Poisson
problems in two and three dimensions. The scheme is based on reformulating the
original problem posed over $\Omega$ on the extruded domain
$\mathcal{C}=\Omega\times[0,\infty)$ following Caffarelli and Silvestre (2007).
The resulting degenerate elliptic integer order PDE is then approximated using
a hybrid FEM-spectral scheme. Finite elements are used in the direction
parallel to the problem domain $\Omega$, and an appropriate spectral method is
used in the extruded direction. The spectral part of the scheme requires that
we approximate the true eigenvalues of the integer order Laplacian over
$\Omega$. We derive an a priori error estimate which takes account of the error
arising from using an approximation in place of the true eigenvalues. We
further present a strategy for choosing approximations of the eigenvalues based
on Weyl's law and finite element discretizations of the eigenvalue problem. The
system of linear algebraic equations arising from the hybrid FEM-spectral
scheme is decomposed into blocks which can be solved effectively using standard
iterative solvers such as multigrid and conjugate gradient. Numerical examples
in two and three dimensions show that the approach is quasi-optimal in terms of
complexity.
