  An additive Runge-Kutta method is used for the time stepping, which
integrates the linear stiff terms by an explicit singly diagonally implicit
Runge-Kutta (ESDIRK) method and the nonlinear terms by an explicit Runge-Kutta
(ERK) method. In each time step, the implicit solve is performed by the
recently developed Hierarchical Poincar\'e-Steklov (HPS) method. This is a fast
direct solver for elliptic equations that decomposes the space domain into a
hierarchical tree of subdomains and builds spectral collocation solvers locally
on the subdomains. These ideas are naturally combined in the presented method
since the singly diagonal coefficient in ESDIRK and a fixed time-step ensures
that the coefficient matrix in the implicit solve of HPS remains the same for
all time stages. This means that the precomputed inverse can be efficiently
reused, leading to a scheme with complexity (in two dimensions)
$\mathcal{O}(N^{1.5})$ for the precomputation where the solution operator to
the elliptic problems is built, and then $\mathcal{O}(N \log N)$ for the solve
in each time step. The stability of the method is proved for first order in
time and any order in space, and numerical evidence substantiates a claim of
stability for a much broader class of time discretization methods. Numerical
experiments supporting the accuracy of efficiency of the method in one and two
dimensions are presented.
