  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. We introduce a new method to simulate
collisional stellar systems, and validate it by comparison with direct $N$-body
codes on small-$N$ simulations. 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. 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. 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 Cen and M54 and even the nuclear
star cluster, beyond the limits of current direct N-body codes in terms of the
number of particles.
