The stochastic model of the growth of an ensemble of GaN nanowires (NW) in a plasma-assisted molecular beam epitaxy suggested in our recent paper (Sabelfeld and Kablukova, 2016) is here further developed to include coalescence caused by bundling of nanowires. Moreover, in the extended model the Ga and N atoms are simulated separately to mimic nucleation of stable GaN islands which then start to grow in one direction with a fixed diameter. The bundling phenomenon is driven by the gain of surface energy at the expense of the elastic energy of bending and becomes energetically favorable once the nanowires exceed a certain critical length as found and described in our paper (Kaganer et al., 2016). The simulation model is based on a probabilistic description of surface diffusion, and takes into account the shading, multiple rescattering of atoms, and their survival probability. The model is implemented in a form of a direct simulation Monte Carlo algorithm. We present a comprehensive analysis of the kinetics of NW growth from an initial height distribution around tens of nanometers to NW heights up to several thousands nanometers which corresponds to the growth time in physical growth experiments about of 3–4 h. We compare the simulation results with our recently developed model without bundling and experimental results as well which show a good agreement. The series of simulations have basically confirmed the remarkable time evolution of the nanowire height distribution, namely, that under some conditions, now described more precisely, the NW height distribution is self-preserving and may become even narrower independent of the form of the initial NW height distribution. This phenomenon is explained by the multiple rescattering of atoms between the NW surfaces.