  We present a simple yet accurate method to compute the adjoint double layer
potential, which is used to solve the Neumann boundary value problem for
Laplace's equation in three dimensions. An expansion in curvilinear coordinates
leads us to modify the expression for the adjoint double layer so that the
singularity is reduced when evaluating the integral on the surface. We then
regularize the Green's function, with a radial parameter $\delta$. We show that
a natural regularization has error $O(\delta^3)$, and a simple modification
improves the error to $O(\delta^5)$. The integral is evaluated numerically
without the need of special coordinates. We use this treatment of the adjoint
double layer to solve the classical integral equation for the interior Neumann
problem and evaluate the solution on the boundary. Choosing $\delta =
ch^{4/5}$, we find about $O(h^4)$ convergence in our examples, where $h$ is the
spacing in a background grid.
