  We introduce a novel spectral element method based on the ultraspherical
spectral method and the hierarchical Poincar\'{e}-Steklov scheme for solving
second-order linear partial differential equations on polygonal domains with
unstructured quadrilateral or triangular meshes. Properties of the
ultraspherical spectral method lead to almost banded linear systems, allowing
the element method to be competitive in the high-polynomial regime ($p > 5$).
The hierarchical Poincar\'{e}-Steklov scheme enables precomputed solution
operators to be reused, allowing for fast elliptic solves in implicit and
semi-implicit time-steppers. The resulting spectral element method achieves an
overall computational complexity of $\mathcal{O}(p^4/h^3)$ for mesh size $h$
and polynomial order $p$, enabling $hp$-adaptivity to be efficiently performed.
We develop an open-source software system, ultraSEM, for flexible,
user-friendly spectral element computations in MATLAB.
