  Second-order gravitational self-force theory has recently led to the
breakthrough calculation of "first post-adiabatic" (1PA) compact-binary
waveforms~[Phys. Rev. Lett. 130, 241402 (2023)]. The computations underlying
those waveforms depend on a method of solving the perturbative second-order
Einstein equation in the Fourier domain. In this paper we present that method,
which involves dividing the domain into several regions. Different regions
utilize different time slicings and allow for the use of "punctures" to tame
sources and enforce physical boundary conditions. We demonstrate the method for
Lorenz-gauge and Teukolsky equations in the relatively simple case of
calculating parametric derivatives ("slow time derivatives") of first-order
fields, which are an essential input at second order.
