We present a new parallel algorithm for solving the Poisson equation in the context of non-stationary stellar dynamics problems, e.g. rotating galaxies or circumstellar disks. This allows us to conduct numerical experiments on a mesh with 10–100 billion of nodes and to use more than 10 thousand of processors. This algorithm is based on a finite difference method, 3D domain decomposition, partial pre-computation of Green's function for subdomains boundaries and parallel implementation of the James method. All the computations have been conducted in the Siberian Supercomputer Center (Novosibirsk) and the Joint Supercomputer Center (Moscow).