  We present a novel implementation of an iterative solver for the solution of
the Poisson equation in the PLUTO code for astrophysical fluid dynamics. Our
solver relies on a relaxation method in which convergence is sought as the
steady-state solution of a parabolic equation, whose time-discretization is
governed by the \textit{Runge-Kutta-Legendre} (RKL) method. Our findings
indicate that the RKL-based Poisson solver, which is both fully parallel and
rapidly convergent, has the potential to serve as a practical alternative to
conventional iterative solvers such as the \textit{Gauss-Seidel} (GS) and
\textit{successive over-relaxation} (SOR) methods. Additionally, it can
mitigate some of the drawbacks of these traditional techniques. We incorporate
our algorithm into a multigrid solver to provide a simple and efficient gravity
solver that can be used to obtain the gravitational potentials in
self-gravitational hydrodynamics. We test our implementation against a broad
range of standard self-gravitating astrophysical problems designed to examine
different aspects of the code. We demonstrate that the results match
excellently with the analytical predictions (when available), and the findings
of similar previous studies.
