  We present a systematic algorithm for the perturbative computation of soft
functions that are defined in terms of two light-like Wilson lines. Our method
is based on a universal parametrisation of the phase-space integrals, which we
use to isolate the singularities in Laplace space. The observable-dependent
integrations can then be performed numerically, and they are implemented in the
new, publicly available package SoftSERVE that we use to derive all of our
numerical results. Our algorithm applies to both SCET-1 and SCET-2 soft
functions, and in the current version it can be used to compute two out of
three NNLO colour structures associated with the so-called correlated-emission
contribution. We confirm existing two-loop results for about a dozen $e^+e^-$
and hadron-collider soft functions, and we obtain new predictions for the
C-parameter as well as thrust-axis and broadening-axis angularities.
