I would like to find a solution to Poisson's equation (in 2D), using the separation of variables. The problem is the following (assuming that $x$ and $y$ are reals)
$u_{xx}+u_{yy}=f(x,y)$ in the rectangle $\Omega=\{0<x<a, 0<y<b\}$
$\frac{\partial u}{\partial \mathbf{n}}=0$ on $\partial \Omega$
$f(x,y) = \left\{ \begin{array}{ll} {-}K & : c<x<d, e<y<f\\ 0 & : elsewhere \end{array} \right.$
where $0<c<d<a$ and $0<e<f<b$.
I am able to solve this problem (with the separation of variables) if one side of the rectangle $\Omega$ is a Dirichlet boundary condition ($u=0$), but not if all 4 sides are Neumann.
In that case (if 3 sides are Neumann and 1 is Dirichlet), the solution is of the form
$u=\sum_{h=1}^{\infty}(k_1 \cos mx + k_2 \sin mx)(k_3 \cosh my + k_4 \sinh my) + k_5 x + k_6 y + k_7$
with $m=\frac{h 2\pi}{2a}$ and the source term is expressed as a Fourier series $K(x)=\sum_{h=1}^\infty K_m \cos mx$.
Perhaps this is where I am wrong when solving this problem with 4 Neumann boundary conditions. I use the same general solution for $u$ (for the complimentary function) and the same particular integral (found using the Fourier expansion for $K$), but try to calculate the new coefficients $k_j$ for the new boundary conditions. Is this how I should proceed? Or would the general solution and Fourier expansion be different for that case?
The pure Neumann problem (when you have the Neumann boundary condition on all parts of the boundary) has two features that distinguish it from the Dirichlet or mixed problems:
To see where $\int_\Omega f=0$ comes from, recall that the Laplacian is the divergence of the gradient field $\nabla u$. By the divergence (Gauss) theorem, $\int_\Omega \Delta u$ is the flux of $\nabla u$ across $\partial \Omega$: $$\int_\Omega \Delta u = \int_{\partial \Omega} \frac{\partial u}{\partial \mathbf n}\tag1$$ (using the outward normal vector). The Neumann condition requires the right-hand side of (1) to be zero.
In your situation, $\int_\Omega f=-K(d-c)(f-e)\ne 0$, hence there is no solution.