  We present an accurate and efficient method to calculate the gravitational
potential of an isolated system in three-dimensional Cartesian and cylindrical
coordinates subject to vacuum (open) boundary conditions. Our method consists
of two parts: an interior solver and a boundary solver. The interior solver
adopts an eigenfunction expansion method together with a tridiagonal matrix
solver to solve the Poisson equation subject to the zero boundary condition.
The boundary solver employs James's method to calculate the boundary potential
due to the screening charges required to keep the zero boundary condition for
the interior solver. A full computation of gravitational potential requires
running the interior solver twice and the boundary solver once. We develop a
method to compute the discrete Green's function in cylindrical coordinates,
which is an integral part of the James algorithm to maintain second-order
accuracy. We implement our method in the {\tt Athena++} magnetohydrodynamics
code, and perform various tests to check that our solver is second-order
accurate and exhibits good parallel performance.
