  We present a new algorithm for numerical computation of large eigenvalues and
associated eigenfunctions of the Dirichlet Laplacian in a smooth, star-shaped
domain in $\mathbb{R}^d$, $d\ge 2$. Conventional boundary-based methods require
a root-search in eigenfrequency $k$, hence take $O(N^3)$ effort per eigenpair
found, using dense linear algebra, where $N=O(k^{d-1})$ is the number of
unknowns required to discretize the boundary. Our method is O(N) faster,
achieved by linearizing with respect to $k$ the spectrum of a weighted interior
Neumann-to-Dirichlet (NtD) operator for the Helmholtz equation. Approximations
$\hat{k}_j$ to the square-roots $k_j$ of all O(N) eigenvalues lying in $[k -
\epsilon, k]$, where $\epsilon=O(1)$, are found with $O(N^3)$ effort. We prove
an error estimate $$ |\hat k_j - k_j| \leq C \Big(\frac{\epsilon^2}{k} +
\epsilon^3 \Big), $$ with $C$ independent of $k$. We present a higher-order
variant with eigenvalue error scaling empirically as $O(\epsilon^5)$ and
eigenfunction error as $O(\epsilon^3)$, the former improving upon the 'scaling
method' of Vergini--Saraceno. For planar domains ($d=2$), with an assumption of
absence of spectral concentration, we also prove rigorous error bounds that are
close to those numerically observed. For $d=2$ we compute robustly the spectrum
of the NtD operator via potential theory, Nystr\"{o}m discretization, and the
Cayley transform. At high frequencies (400 wavelengths across), with
eigenfrequency relative error $10^{-10}$, we show that the method is $10^3$
times faster than standard ones based upon a root-search.
