  Precision studies of the Higgs boson at future $e^+e^-$ colliders can help to
shed light on fundamental questions related to electroweak symmetry breaking,
baryogenesis, the hierarchy problem, and dark matter. The main production
process, $e^+e^- \to HZ$, will need to be controlled with sub-percent
precision, which requires the inclusion of next-to-next-to-leading order (NNLO)
electroweak corrections. The most challenging class of diagrams are planar and
non-planar double-box topologies with multiple massive propagators in the
loops. This article proposes a technique for computing these diagrams
numerically, by transforming one of the sub-loops through the use of Feynman
parameters and a dispersion relation, while standard one-loop formulae can be
used for the other sub-loop. This approach can be extended to deal with tensor
integrals. The resulting numerical integrals can be evaluated in minutes on a
single CPU core, to achieve about 0.1% relative precision.
