  Exceptionally elegant formulae exist for the fractional Laplacian operator
applied to weighted classical orthogonal polynomials. We utilize these results
to construct a solver, based on frame properties, for equations involving the
fractional Laplacian of any power, $s \in (0,1)$, on an unbounded domain in one
or two dimensions. The numerical method represents solutions in an expansion of
weighted classical orthogonal polynomials as well as their unweighted
counterparts with a specific extension to $\mathbb{R}^d$, $d \in \{1,2\}$. We
examine the frame properties of this family of functions for the solution
expansion and, under standard frame conditions, derive an a priori estimate for
the stationary equation. Moreover, we prove one achieves the expected order of
convergence when considering an implicit Euler discretization in time for the
fractional heat equation. We apply our solver to numerous examples including
the fractional heat equation (utilizing up to a $6^\text{th}$-order
Runge--Kutta time discretization), a fractional heat equation with a
time-dependent exponent $s(t)$, and a two-dimensional problem, observing
spectral convergence in the spatial dimension for sufficiently smooth data.
