  The Fokker-Planck (FP) model is one of the commonly used methods for studies
of the dynamical evolution of dense spherical stellar systems such as globular
clusters and galactic nuclei. The FP model is numerically stable in most cases,
but we find that it encounters numerical difficulties rather often when the
effects of tidal shocks are included in two-dimensional (energy and angular
momentum space) version of the FP model or when the initial condition is
extreme (e.g., a very large cluster mass and a small cluster radius). To avoid
such a problem, we have developed a new integration scheme for a
two-dimensional FP equation by adopting an Alternating Direction Implicit (ADI)
method given in the Douglas-Rachford split form. We find that our ADI method
reduces the computing time by a factor of ~2 compared to the fully implicit
method, and resolves problems of numerical instability.
