  In this paper, we introduce a novel method for deriving higher order
corrections to the mean-field description of the dynamics of interacting
bosons. More precisely, we consider the dynamics of $N$ $d$-dimensional bosons
for large $N$. The bosons initially form a Bose-Einstein condensate and
interact with each other via a pair potential of the form
$(N-1)^{-1}N^{d\beta}v(N^\beta\cdot)$ for $\beta\in[0,\frac{1}{4d})$. We derive
a sequence of $N$-body functions which approximate the true many-body dynamics
in $L^2(\mathbb{R}^{d N})$-norm to arbitrary precision in powers of $N^{-1}$.
The approximating functions are constructed as Duhamel expansions of finite
order in terms of the first quantised analogue of a Bogoliubov time evolution.
