How to derive this? (Green's Second Identity, Poisson Boltzmann Equation)

I'm trying to derive the integral form of $\phi_p^{ext}$ and $\phi_p^{int}$ as in this paper: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2693949/.

Specifically, I'm trying to go from (2) to (3), (4) in their equation numberings in section 2.

They apply Green's second identity to the linearized Poisson Boltzmann equation to get the integral form of $\phi_p^{ext}$ and $\phi_p^{int}$.

This is the linearized PB equation (2): $$-\nabla(\epsilon\nabla\phi)+\kappa^2\phi = \sum_{i=1}^{M}q_i\delta(r-r_i)$$

I'd appreciate any help! There are two things that aren't immediately obvious from the equations, that on the interface between the interior ($int$) and exterior ($ext$), $\phi_p^{ext}=\phi_p^{int}$ and $\frac{\partial\phi^{int}}{\partial\mathbf{n}} = \frac{\partial\phi^{ext}}{\partial\mathbf{n}}$. These are written up in the paragraph between (4) and (5) in the linked paper, written here for convenience.

The Green function $G_{pt}$ satisfies $$-\nabla(\epsilon\nabla G_{pt})+\kappa^2 G_{pt}=\delta(t-p),$$ where $\nabla$ describes derivatives with respect to $t$. Physically, that is, $G_{pt}$ is the field at the point $t$ due to a point source at the point $p$, under some (unspecified) set of boundary conditions. Multiply the left side of the linearized PDE, as a function of $t$, by $G_{pt}$, and integrate over $\Omega$. Moving the derivatives over to act on the Green function via two integrations by parts gives the surface terms shown, and a remaining volume integral of $\delta(t-p)\, \phi^{\mathrm{int}}(t)$.
The game in the paper seems to be that, where the potential is large and quickly varying, i.e., near the cluster of sources that make up the molecule, the sine term in the original, non-linear PDE is negligible, and one has effectively a Poisson equation. (Or maybe the point is that the fluid has different properties on either side of a physical boundary. Guess who's not a chemist.) Farther away, where the potential is small, one can linearize and get a (modified) Helmholtz equation. Imposing $\phi_t^{\mathrm{int}}=\phi_t^{\mathrm{ext}}$ ties the solutions of those two different PDEs together at the boundary separating "near" from "far."
Eqs. (5) and (6) arise by taking the limit as $p$ approaches the boundary from the inside and outside, respectively. (The principal value gets around the singularity in $G_{pt}$, say, at $p = t$ on the surface.) The boundary values of the potential and electric field need to be related by the two equations in the text just before eq. (5). So, one has two integral equations giving the boundary potential at one point $p$ on the surface, now called $f_p$, in terms of its values $f_t$ everywhere on the surface and the values $h_t$ of the (exterior) normal derivative. Solving these coupled integral equations, which one can hope to do numerically by formulating them as a discretized linear algebra problem, gives a consistent set of boundary data. These boundary data then give the solutions everywhere via eqs. (3) and (4).
So, it shouldn't be obvious that $\phi_p^{\mathrm{int}} = \phi_p^{\mathrm{ext}}$ on the surface. Demanding that equality, and the other relation between the normal derivatives, * is equivalent to solving the PDE*, and must be done numerically.