We prove that the flow of the discrete nonlinear Schrödinger equation (DNLS) is the mean field limit of the quantum dynamics of the Bose–Hubbard model for N interacting particles. In particular, we show that the Wick symbol of the annihilation operators evolved in the Heisenberg picture converges, as N becomes large, to the solution of the DNLS. A quantitative $$L^p$$
L
p
-estimate, for any $$p \ge 1$$
p
≥
1
, is obtained with a linear dependence on time due to a Gaussian measure on initial data coherent states.