In this paper, we present our conclusions from the numerical study of the collapse of a destabilized collisionless stellar system. We use both direct integration of the Vlasov-Poisson equations and an N-body tree code to obtain our results, which are mutually confirmed. We find that spherical and moderately nonspherical collapse configurations evolve to new equilibrium configurations in which the velocity distribution approaches a Gaussian form, at least in the central regions. The evolution to this state has long been an open question, and in this work we are able to clarify the process responsible and to support predictions made from statistical considerations (Lynden-Bell 1967; Nakamura 2000). The simulations of merging N-body systems show a transition to a Gaussian velocity distribution that is increasingly suppressed as the initial separation of centres is increased. Possible reasons for this are discussed.