  We consider a system of N interacting bosons in the mean-field scaling regime
and construct corrections to the Bogoliubov dynamics that approximate the true
N-body dynamics in norm to arbitrary precision. The N-independent corrections
are given in terms of the solutions of the Bogoliubov and Hartree equations and
satisfy a generalized form of Wick's theorem. We determine the n-point
correlation functions of the excitations around the condensate, as well as the
reduced densities of the N-body system, to arbitrary accuracy, given only the
knowledge of the two-point functions of a quasi-free state and the solution of
the Hartree equation. In this way, the complex problem of computing all n-point
correlation functions for an interacting N-body system is essentially reduced
to the problem of solving the Hartree equation and the PDEs for the Bogoliubov
two-point functions.
