Derivation of a series expansion used to solve Laplace's equation The paper Mathematics of the Faraday Cage contains an interesting discussion of the Faraday cage effect. I am fine with the physics, and they appear to have done a good job with the numerics. I am however baffled by the form of the series expansion they have used to set up the numerical problem. I will present the form of the Laplace equation as it appears in the paper, and then present the series expansion which confuses me.
The problem is as follows: consider the complex plane with $n$ disks of radius $r$ centered at the $n$ roots of unity.

The authors seek a real function $\phi(z)$ that satisfies the Laplace equation  
$\nabla^2 \phi = 0$  
in the region in the picture above exterior to the unit circle $\Gamma$, and the boundary condition $\phi=\phi_0$ on the disks ($\phi_0$ is unknown and part of the solution). 
They also require $\phi(z)$ to have the asymptotic behavior  
$\phi(z) = \log|z-z_s| + O(1)$ as $z \rightarrow z_s$,  
$\phi(z) = \log|z| + o(1)$ as $z \rightarrow \infty$,
where $z_s$ is a fixed point external to the 'cage' that will serve as the location of a test charge.
The series expansion they use is: 
$\phi(z) = \log|z -z_s| + \sum_{j=1}^n \left(d_j\log|z-c_j| + \mathrm{Re}\left[\sum_{k=1}^N({a_{jk}}-ib_{jk})(z-c_j)^{-k} \right]\right)$.
In this expansion $N$ is free, and determines the accuracy of the numerical solution. The sets $\{d_j
\}$, $\{a_{jk}\}$ and $\{b_{jk}\}$ are real constants to be determined, and the $\{c_j\}$ are the centers of the $n$ disks
Can anyone help me understand how this expansion was derived? Or at least why it is valid?
Chapman, S.Jonathan; Hewett, David P.; Trefethen, Lloyd N., Mathematics of the Faraday cage, SIAM Rev. 57, No. 3, 398-417 (2015). ZBL1339.31001.. 
 A: I managed to figure this out by rereading Appendix A in the paper I linked to (the comment of @QtizedD was also useful, thanks)
Firstly, it is obvious that the series expansion satisfies $\nabla^2 \phi = 0$: the simplest way to see this is probably to use the Wirtinger derivative representation of the Laplacian, $\nabla^2 u = 4 \frac{\partial^2 u}{\partial z \partial \bar{z}}$. The specific choice of the matrix coefficients $a_{ij}$, $b_{ij}$ and the $d_j$ are necessary to satisfy the additional requirement that there be constant potential on the wires, $\phi = \phi_0$.
The first log term is clearly the 2D Green’s function for the source charge. The other log terms are the Green’s functions for the disks on the unit circle, with the restraint that the total charge on all the disks is zero ($\sum d_j = 0$), and with the constant potential. 
The step I was missing was that, in addition to choosing $N$ to give us the desired accuracy in the numerical solution, the disks are themselves discretised, yielding an overdetermined linear system, that can be numerically approximated by the method of least squares.
As a bonus, here is a picture I generated with Python for the above, for 10 disks with a radius of 0.1, and the source charge located at -2 (the black contours are the potential $\phi(z)$, and the red dot is the source charge):

