How is Neumann's green function found for a particular surface? Green's Functions can be used to solve boundary value problems. In particular, there is a textbook standard recipe for dealing with von Neumann boundary conditions.
But how can I find a particular Green's function?
Example:
A cuboid width=a, length=b and height=c, defining a simple surface.
According to textbooks (e.g. Jackson, Classical Electrodynamics) Neumann's Green function $G_n(x,x')$ is defined by
$\Delta' G_n(x,x') = -4 \pi \delta(x-x')$
and
$\partial G / \partial n' = -4 \pi / S = const.$   
for x' on the surface
Here S is the total surface of the cuboid.
I understand the idea behind completely, but have no Idea, how $G_n(x,x')$ is to be constructed for the particular geometry of the example.
 A: Hint:
Try the method of image charges, augmenting the original box to a rectangular lattice (one with reflection symmetries, so the field vanishes across each face.)

More explicit hint:
Intuitively, the idea is to place charges around the boundary $S$ of the cuboid in order to impose von Neumann b.c. (vanishing electric field on $S$.) In general this is difficult, but in this particular case there are essentially two symmetries we can use: reflection and translation.
First, if a distribution of charges has a reflection symmetry about some plane $\hat n$ (i.e. if $\rho(x)=\rho(x-2(x\cdot\hat n)\hat n)$), then the corresponding electric field vanishes on the plane of symmetry (the set of points with $x\cdot \hat n=0$.) If there are multiple planes of reflection symmetry, then the electric field vanishes along each. 
Based on this, it might seem reasonable to define image charges in the following pattern:

Here, the original domain contains a point source (blue circle), and image charges are placed to create a lattice of planes of reflection symmetry. There is one problem with this arrangement, however: in the language of physics, Gauss' law tells us that $\oint_S \hat n(A)\cdot E(A)dA=\frac{1}{\epsilon_0}Q_{\text{enclosed}}\neq 0$, but $E(A)$ must vanish when von Neumann boundary conditions are imposed. To cancel this, we fill space uniformly with a negative background charge to cancel the net positive charge in each block. You can then derive an exact series expression for the Green function (in 'real' or 'Fourier' space) that can be approximated to arbitrary accuracy.
