  We consider the numerical solution of the real time equilibrium Dyson
equation, which is used in calculations of the dynamical properties of quantum
many-body systems. We show that this equation can be written as a system of
coupled, nonlinear, convolutional Volterra integro-differential equations, for
which the kernel depends self-consistently on the solution. As is typical in
the numerical solution of Volterra-type equations, the computational bottleneck
is the quadratic-scaling cost of history integration. However, the structure of
the nonlinear Volterra integral operator precludes the use of standard fast
algorithms. We propose a quasilinear-scaling FFT-based algorithm which respects
the structure of the nonlinear integral operator. The resulting method can
reach large propagation times, and is thus well-suited to explore quantum
many-body phenomena at low energy scales. We demonstrate the solver with two
standard model systems: the Bethe graph, and the Sachdev-Ye-Kitaev model.
