  We present a simple numerical scheme for perturbation theory (PT)
calculations of large-scale structure. Solving the evolution equations for
perturbations numerically, we construct the PT kernels as building blocks of
statistical calculations, from which the power spectrum and/or correlation
function can be systematically computed. The scheme is especially applicable to
the generalized structure formation including modified gravity, in which the
analytic construction of PT kernels is intractable. As an illustration, we show
several examples for power spectrum calculations in $f(R)$ gravity and
$\Lambda$CDM models.
