In this talk I will introduce some of the new features of pySecDec, an automated tool for the numerical computation of Feynman integrals. In particular, I will discuss the automatic application of the method of "expansion by regions" to Feynman integrals and a new, publicly available, feature for numerically evaluating amplitudes. As an example use case, I will present the two-loop virtual corrections to gg→ZH involving a closed top-quark loop. These amplitudes provide the missing ingredient required to compute the NLO QCD corrections to Higgs boson production in association with a Z-boson in gluon fusion. I will describe the techniques used to obtain and evaluate these amplitudes and present comparisons to various approximations in the literature.