I have been wrestling with the following problem for several days to get an analytical solution. Any suggestions or guideline will be greatly appreciated. The PDE is being solved over a cube, $ -a\le x \le a$; $-b\le y \le b$; $0 \le z \le L; $
$ \nabla ^2 u=f(x,y) $
$ {\partial u \over \partial x}=0 $ $\;\;\;\;$at $ x=-a$ and $x=a$;
$ {\partial u \over \partial y}=0 $ $\;\;\;\;$at $ y=-b$ and $y=b$;
$ {\partial u \over \partial z}=Cx $ $\;$at $ z=0$ where $C$ is a constant;
$ {\partial u \over \partial z}=0 $ $\;\;\;\;$at $ z=L$
Hint: Take $\partial/\partial z$ of both sides to obtain
$$ \nabla^2u_z = 0 $$
Perform separation of variables by setting $u_z(x,y,z) = X(x)Y(y)Z(z)$ to get
$$ \frac{X''}{X} + \frac{Y''}{Y} = -\frac{Z''}{Z} = -\lambda^2 $$
The third boundary condition implies that the $X$ part must be odd and the $Y$ part must be even, so we have the following solution forms
$$ X_n = \sin\left[\left(n+\frac{1}{2}\right)\frac{\pi x}{a}\right] $$
$$ Y_m = \cos\left( m\frac{\pi y}{b} \right) $$
$$ Z_{nm} = \frac{\sinh\left[\lambda_{nm}(L-z)\right]}{\sinh(\lambda_{nm} L)} $$
where
$$ \lambda_{nm}^2 = \left[\left(n+\frac{1}{2}\right)\frac{\pi}{a}\right]^2 + \left( m\frac{\pi}{b} \right)^2 $$