  We report a numerical calculation of the two-photon absorption coefficient of
electrons in a binding potential using the real-time real-space higher-order
difference method. By introducing random vector averaging for the intermediate
state, the task of evaluating the two-dimensional time integral is reduced to
calculating two one-dimensional integrals. This allows the reduction of the
computation load down to the same order as that for the linear response
function. The relative advantage of the method compared to the straightforward
multi-dimensional time integration is greater for the calculation of non-linear
response functions of higher order at higher energy resolution.
