Non-linear effects are crucial for computing the cosmological matter power spectrum to theaccuracy required by future generation surveys. Here, a new approach is presented,in which the power spectrum, the bispectrum and higher order correlations areobtained—at any redshift and for any momentum scale—by integrating a systemof differential equations. The method is similar to that of the familiar BBGKY(Bogoliubov–Born–Green–Kirkwood–Yvon) hierarchy. Truncating at the level ofthe trispectrum, the solution of the equations corresponds to the summation ofan infinite class of perturbative corrections. Compared to other resummationframeworks, the scheme discussed here is particularly suited to cosmologies other thanΛCDM (CDM: cold dark matter), such as those based on modifications of gravity and those containingmassive neutrinos. As a first application, we compute the baryonic acoustic oscillation featureof the power spectrum, and compare the results with perturbation theory, the halo model, andN-body simulations. The density–velocity and velocity–velocity power spectra are alsocomputed, revealing that they are much less contaminated by non-linearities than thedensity–density one. The approach can be seen as a particular formulation of therenormalization group, in which time is the flow parameter.