  We provide a novel approach to the numerical solution of the family of
nonlocal elliptic equations $(-\Delta)^su=f$ in $\Omega$, subject to some
homogeneous boundary conditions $\mathcal{B}(u)=0$ on $\partial \Omega$, where
$s\in(0,1)$, $\Omega\subset \mathbb{R}^n$ is a bounded domain, and
$(-\Delta)^s$ is the spectral fractional Laplacian associated to $\mathcal{B}$
on $\partial \Omega$. We use the solution representation $(-\Delta)^{-s}f$
together with its singular integral expression given by the method of
semigroups. By combining finite element discretizations for the heat semigroup
with monotone quadratures for the singular integral we obtain accurate
numerical solutions. Roughly speaking, given a datum $f$ in a suitable
fractional Sobolev space of order $r\geq 0$ and the discretization parameter
$h>0$, our numerical scheme converges as $O(h^{r+2s})$, providing super
quadratic convergence rates up to $O(h^4)$ for sufficiently regular data, or
simply $O(h^{2s})$ for merely $f\in L^2(\Omega)$. We also extend the proposed
framework to the case of nonhomogeneous boundary conditions and support our
results with some illustrative numerical tests.
