  The numerical simulation of weakly nonlinear ultrasound is important in
treatment planning for focused ultrasound (FUS) therapies. However, the large
domain sizes and generation of higher harmonics at the focus make these
problems extremely computationally demanding. Numerical methods typically
employ a uniform mesh fine enough to resolve the highest harmonic present in
the problem, leading to a very large number of degrees of freedom. This paper
proposes a more efficient strategy in which each harmonic is approximated on a
separate mesh, the size of which is proportional to the wavelength of the
harmonic. The increase in resolution required to resolve a smaller wavelength
is balanced by a reduction in the domain size. This nested meshing is feasible
owing to the increasingly localised nature of higher harmonics near the focus.
  Numerical experiments are performed for FUS transducers in homogeneous media
in order to determine the size of the meshes required to accurately represent
the harmonics. In particular, a fast \emph{volume potential} approach is
proposed and employed to perform convergence experiments as the computation
domain size is modified. This approach allows each harmonic to be computed via
the evaluation of an integral over the domain. Discretising this integral using
the midpoint rule allows the computations to be performed rapidly with the FFT.
It is shown that at least an order of magnitude reduction in memory consumption
and computation time can be achieved with nested meshing. Finally, it is
demonstrated how to generalise this approach to inhomogeneous propagation
domains.
