Context : Stellar systems are broadly divided into collisional and non-collisional .
The latter are large- N systems with long relaxation timescales and can be simulated disregarding two-body interactions , while either computationally expensive direct N -body simulations or approximate schemes are required to properly model the former .
Large globular clusters and nuclear star clusters , with relaxation timescales of the order of a Hubble time , are small enough to display some collisional behaviour and big enough to be impossible to simulate with direct N -body codes and current hardware .

Aims : We introduce a new method to simulate collisional stellar systems , and validate it by comparison with direct N -body codes on small-N simulations .

Methods : The Multi-Particle collision for Dense stellar systems Code ( mpcdss ) is a new code for evolving stellar systems with the Multi-Particle Collision method .
Such method amounts to a stochastic collision rule that allows to conserve exactly the energy and momentum over a cluster of particles experiencing the collision .
The code complexity scales with N \log N in the number of particles .
Unlike Monte-Carlo codes , mpcdss can easily model asymmetric , non-homogeneous , unrelaxed and rotating systems , while allowing us to follow the orbits of individual stars .

Results : We evolve small ( N = 3.2 \times 10 ^ { 4 } ) star clusters with mpcdss and with the direct-summation code nbody6 , finding a similar evolution of key indicators .
We then simulate different initial conditions in the 10 ^ { 4 } -10 ^ { 6 } star range .

Conclusions : mpcdss bridges the gap between small , collisional systems that can be simulated with direct N -body codes and large noncollisional systems . mpcdss in principle allows us to simulate globular clusters such as Omega Centauri and M54 and even the nuclear star cluster , beyond the limits of current direct N -body codes in terms of the number of particles .