  We have developed a novel Monte Carlo method for simulating the dynamical
evolution of stellar systems in arbitrary geometry. The orbits of stars are
followed in a smooth potential represented by a basis-set expansion and
perturbed after each timestep using local velocity diffusion coefficients from
the standard two-body relaxation theory. The potential and diffusion
coefficients are updated after an interval of time that is a small fraction of
the relaxation time, but may be longer than the dynamical time. Thus our
approach is a bridge between the Spitzer's formulation of the Monte Carlo
method and the temporally smoothed self-consistent field method. The primary
advantages are the ability to follow the secular evolution of shape of the
stellar system, and the possibility of scaling the amount of two-body
relaxation to the necessary value, unrelated to the actual number of particles
in the simulation. Possible future applications of this approach in galaxy
dynamics include the problem of consumption of stars by a massive black hole in
a non-spherical galactic nucleus, evolution of binary supermassive black holes,
and the influence of chaos on the shape of galaxies, while for globular
clusters it may be used for studying the influence of rotation.
