  We present an approach for simulating the collisional evolution of spherical
isotropic stellar systems based on the one-dimensional Fokker-Planck equation.
A novel aspect is that we use the phase volume as the argument of the
distribution function, instead of the traditionally used energy, which
facilitates the solution. The publicly available code, PhaseFlow, implements a
high-accuracy finite-element method for the Fokker-Planck equation, and can
handle multiple-component systems, optionally with the central black hole and
taking into account loss-cone effects and star formation. We discuss the energy
balance in the general setting, and in application to the Bahcall-Wolf cusp
around a central black hole, for which we derive a perturbative solution. We
stress that the cusp is not a steady-state structure, but rather evolves in
amplitude while retaining an approximately $\rho\propto r^{-7/4}$ density
profile. Finally, we apply the method to the nuclear star cluster of the Milky
Way, and illustrate a possible evolutionary scenario in which a two-component
system of lighter main-sequence stars and stellar-mass black holes develops a
Bahcall-Wolf cusp in the heavier component and a weaker $\rho\propto r^{-3/2}$
cusp in the lighter, visible component, over the period of several Gyr. The
present-day density profile is consistent with the recently detected mild cusp
inside the central parsec, and is weakly sensitive to initial conditions.
