  A spectral method for solving linear partial differential equations (PDEs)
with variable coefficients and general boundary conditions defined on
rectangular domains is described, based on separable representations of partial
differential operators and the one-dimensional ultraspherical spectral method.
If a partial differential operator is of splitting rank $2$, such as the
operator associated with Poisson or Helmholtz, the corresponding PDE is solved
via a generalized Sylvester matrix equation, and a bivariate polynomial
approximation of the solution of degree $(n_x,n_y)$ is computed in
$\mathcal{O}((n_x n_y)^{3/2})$ operations. Partial differential operators of
splitting rank $\geq 3$ are solved via a linear system involving a block-banded
matrix in $\mathcal{O}(\min(n_x^{3} n_y,n_x n_y^{3}))$ operations. Numerical
examples demonstrate the applicability of our 2D spectral method to a broad
class of PDEs, which includes elliptic and dispersive time-evolution equations.
The resulting PDE solver is written in MATLAB and is publicly available as part
of CHEBFUN. It can resolve solutions requiring over a million degrees of
freedom in under $60$ seconds. An experimental implementation in the Julia
language can currently perform the same solve in $10$ seconds.
