We present a clear N-body realization of the growth of
a Bahcall-Wolf f ~ E1/4, rho ~ r-7/4 density cusp around a massive object ("black hole") at the
center of a stellar system. Our N-body algorithm incorporates a novel implementation of
Mikkola-Aarseth chain regularization to handle close interactions between star and black hole particles.
Forces outside the chain were integrated on a GRAPE-6A/8 special-purpose computer with particle numbers up to N=0.25 x 10^6.
We compare our N-body results with predictions of the isotropic Fokker-Planck equation and verify that the time dependence
of the density (both configuration and phase-space) predicted by the Fokker-Planck equation is well reproduced by the N-body
algorithm, for various choices of N and of the black hole mass. Our results demonstrate the feasibility of direct-force integration
techniques for simulating the evolution of galactic nuclei on relaxation time scales.