  The distributions of the $k$-th largest level at the soft edge scaling limit
of Gaussian ensembles are some of the most important distributions in random
matrix theory, and their numerical evaluation is a subject of great practical
importance. One numerical method for evaluating the distributions uses the fact
that they can be represented as Fredholm determinants involving the so-called
Airy integral operator. When the spectrum of the integral operator is computed
by discretizing it directly, the eigenvalues are known to at most absolute
precision. Remarkably, the Airy integral operator is an example of a so-called
bispectral operator, which admits a commuting differential operator that shares
the same eigenfunctions. In this manuscript, we develop an efficient numerical
algorithm for evaluating the eigendecomposition of the Airy integral operator
to full relative precision, using the eigendecomposition of the commuting
differential operator. This allows us to rapidly evaluate the distributions of
the $k$-th largest level to full relative precision rapidly everywhere except
in the left tail, where they are computed to absolute precision. In addition,
we characterize the eigenfunctions of the Airy integral operator, and describe
their extremal properties in relation to an uncertainty principle involving the
Airy transform. We observe that the Airy integral operator is fairly universal,
and we describe a separate application to Airy beams in optics. Using the
eigenfunctions, we compute a finite-energy Airy beam that is optimal, in the
sense that the beam is both maximally concentrated, and maximally
non-diffracting and self-accelerating.
