ddition, recent observations of young stellar objects indicate that at least in the early stage of evolution, protostellar disks are massive enough to be self-gravitating [@kratter16; @tobin16]. Gravitational instability of such disks may form trailing spirals that can redistribute the mass and angular momentum and induce heat via shocks [@mejia05; @evans15] and may be responsible for the formation of giant planets [@boss07; @zhu12].

To follow evolution of self-gravitating disks, one needs to solve the Poisson equation $$\label{eq:Poisson}
  \nabla^2\Phi = 4\pi G \rho,$$ in cylindrical coordinates $(R,\phi,z)$ subject to a proper boundary condition. In Equation , $\Phi$, $\rho$, and $G$ refer to the gravitational potential, mass density, and gravitational constant, respectively. For an isolated system, $\Phi$ has to satisfy vacuum (or “open”) boundary conditions (i.e., $\Phi$ vanishes at infinite distances), for which the formal solution of Equation is given by $$\label{eq:intPoisson}
  \Phi({\bf x}) = \iiint {\cal G}_\infty({\bf x,x'}) \rho({\bf x'}) \,d^3x',$$ where ${\cal G}_\infty({\bf x,x'})  \equiv - G/ |{\bf x-x'}|$ is the gravitational potential per unit mass due to a point source situated at ${\bf x'}$. Hereafter, we call ${\cal G}_\infty({\bf x,x'})$ the continuous Green’s function (CGF) to distinguish it from the discrete Green’s function (DGF) based on the *discrete* Laplace operator (e.g., @burk97) discussed in Section \[s:dgf\].

In simulating dynamics of geometrically thin disks, it has been customary to assume that the disk density along the vertical direction follows a simple function such as Dirac’s delta function for a razor-thin disk and a Gaussian function for a slightly extended disk (e.g., @kal71 [@miller76; @li09; @wang15]). In this case, the integral along the $z$-direction in Equation can be performed analytically, and finding $\Phi(R, \phi)$ at $z=0$ reduces to numerical evaluation of the remaining two-dimensional (2D) integral in the $R$–$\phi$ plane. For example, @miller7