  The vertical modes of linearized equations of motion are widely used by the
oceanographic community in numerous theoretical and observational contexts.
However, the standard approach for solving the generalized eigenvalue problem
using second-order finite difference matrices produces $O(1)$ errors for all
but the few lowest modes, and increasing resolution quickly becomes too slow as
the computational complexity of eigenvalue algorithms increase as $O(n^3)$.
Existing methods are therefore inadequate for computing a full spectrum of
internal waves, such as needed for initializing a numerical model with a full
internal wave spectrum. Here we show that rewriting the eigenvalue problem in
stretched coordinates and projecting onto Chebyshev polynomials results in
substantially more accurate modes than finite-differencing at a fraction of the
computational cost. We also compute the surface quasigeostrophic modes using
the same methods. All spectral and finite difference algorithms are made
available in a suite of Matlab classes that have been validated against known
analytical solutions in constant and exponential stratification.
