The design of phased array transducers for medical diagnostic ultrasound asks for an understanding of the nonlinear propagation of acoustic wavefields. Most existing numerical models are based on the linearized model equations, but in the recent decades several numerical models have been developed that incorporate weak nonlinear propagation as well. In this study, we present an approach that enables the computation of large-scale, three-dimensional, weakly nonlinear wavefields in the time domain. It is based on the Neumann iterative method, which enables the successive use of the solution of a linear wave problem, where the nonlinearity is treated as a contrast source. The linear wave problem is formulated as a convolution integral with a f...