In galactic nuclei, the gravitational potential is dominated by the central supermassive black hole, so stars follow quasi-Keplerian orbits. These orbits are distorted by gravitational forces from other stars, leading to long-term orbital relaxation. The direct numerical study of these processes is challenging because the fast orbital motion imposed by the central black hole requires very small timesteps. Within the secular approximation of smearing out stars along their underlying Keplerian orbits, a multipole expansion of the pairwise interaction between the stars yields an efficient numerical code to investigate the long-term evolution of their orbital parameters. These new simulations precisely recover the diffusion coefficients of stellar eccentricities obtained through analytical calculations of the secular dynamics. The computational complexity of the present method scales linearly with the total number of stars, so it should prove useful to study the long-term evolution of self-gravitating systems dominated by a central mass.