  We introduce a Fourier-Bessel-based spectral solver for Cauchy problems
featuring Laplacians in polar coordinates under homogeneous Dirichlet boundary
conditions. We use FFTs in the azimuthal direction to isolate angular modes,
then perform discrete Hankel transform (DHT) on each mode along the radial
direction to obtain spectral coefficients. The two transforms are connected via
numerical and cardinal interpolations. We analyze the boundary-dependent error
bound of DHT; the worst case is $\sim N^{-3/2}$, which governs the method, and
the best $\sim e^{-N}$, which then the numerical interpolation governs. The
complexity is $O[N^3]$. Taking advantage of Bessel functions being the
eigenfunctions of the Laplacian operator, we solve linear equations for all
times. For non-linear equations, we use a time-splitting method to integrate
the solutions. We show examples and validate the method on the two-dimensional
wave equation, which is linear, and on two non-linear problems: a
time-dependent Poiseuille flow and the flow of a Bose-Einstein condensate on a
disk.
