Constants in Laplace's equation for a cube I'm working a Laplace's equation $\Delta F=0$  for a cube in Cartesian coordinates $((0,0,0),(a,a,a))$ and after separation I have $$\frac{X''(x)}{X(x)} +\frac{Y''(y)}{Y(y)}+\frac{Z''(z)}{Z(z)} = 0\tag{1}$$
with boundary conditions
$$ x=0 \Rightarrow F=F_0$$
$$ x=a \Rightarrow F=-F_0$$
$$ y=0 \Rightarrow F=0$$
$$ y=a \Rightarrow F=0$$
$$ z=0 \Rightarrow F=0$$
$$ z=a \Rightarrow F=0$$
Thinking about my boundary condition, I figured I'd want the Y and Z terms in $(1)$ to be equal to negative constants $-k^2$ and $-\ell^2$ so I could get a $\sin{\pi y\over a}\sin{\pi z\over a}$ out of them, and the X term to equal a negative constant $-j^2$ as well whence I'd take a $F_0\cos{2\pi x\over a}$. Problem is, if they are all negative constants, they can't add up to zero as in $(1)$! What other options do I have here? I don't want to take a positive $-j^2$ because 

Scale the variable such that $x' = \pi x/a$, $y' = \pi y/a$, $z' = \pi z/a$ first, we the equation still holds
$$
\frac{X''(x')}{X(x')} +\frac{Y''(y')}{Y(y')}+\frac{Z''(z')}{Z(z')} = 0
$$
and we can still get what you mentioned in your question:
$$
\begin{cases}
Y'' = -k^2 Y
\\
Y(0) = Y(\pi) = 0
\end{cases},
\quad
\begin{cases}
Z'' = -l^2 Z
\\
Z(0) = Z(\pi) = 0
\end{cases}.
$$
Slightly abuse of notation, these two lead to 
$$
Y_k(y) = \sin(ky),\quad\text{and}\quad Z_l(z) = \sin(lz).
$$
Now $X$ satisfies:
$$
X'' = (k^2 + l^2)X,
$$
hence let $m = \sqrt{k^2+l^2}$:
$$
X = a_{kl} e^{mx} + b_{kl} e^{-mx},
$$
where the coefficient $a_{kl}, b_{kl}$ depend on $k$ and $l$.
When $x = 0$:
$$
F\big|_{x=0} =F_0  = \sum_{l=1}^{\infty}\sum_{k=1}^{\infty} (a_{kl}+b_{kl}) \sin(ky)\sin(lz)\tag{1}
$$
When $x=\pi$:
$$
F\big|_{x=\pi} =-F_0  = \sum_{l=1}^{\infty}\sum_{k=1}^{\infty} (a_{kl} e^{m\pi} + b_{kl} e^{-m\pi}) \sin(ky)\sin(lz) \tag{2}
$$
For $\{ \sin(ky)\sin(lz)\}_{k,l=1}^{\infty}$ are orthogonal pairwise, multiply (1) by $\sin(ky)\sin(lz)$ and integrate on $(0,\pi)$ for $y$ and $z$:
$$
F_0 \int^{\pi}_0\int^{\pi}_0 \sin(ky)\sin(lz)\,dydz =\int^{\pi}_0\int^{\pi}_0  (a_{kl}+b_{kl}) \sin^2(ky)\sin^2(lz)\,dydz, \tag{$1'$}
$$
The right side only have one term for $k' \neq k$ or $l' \neq l$
$$
\int^{\pi}_0\int^{\pi}_0 \sin(ky)\sin(lz) \sin(k'y)\sin(l'z) \,dydz= 0,
$$ 
this is the orthogonality I mentioned above. Now $(1')$ gives:
$$
F_0\frac{1-(-1)^k}{k}\frac{1-(-1)^l}{l} = (a_{kl}+b_{kl})\frac{\pi^2}{4}.
$$
Doing the same thing for (2) yields:
$$
-F_0\frac{1-(-1)^k}{k}\frac{1-(-1)^l}{l} =  (a_{kl} e^{m\pi} + b_{kl} e^{-m\pi})\frac{\pi^2}{4}.
$$
When $k$ and $l$ are both odd, we have:
$$
a_{kl} = - \frac{16F_0 (e^{-m\pi}+1)}{\pi^2 kl(e^{m\pi} - e^{-m\pi})},\text{ and }
b_{kl} = \frac{16F_0 (e^{m\pi}+1)}{\pi^2 kl(e^{m\pi} - e^{-m\pi})},
$$
When one of $k$ and $l$ is even, we have:
$$
a_{kl} = b_{kl} = 0.
$$
Now the final solution can be obtained just by scaling back by the factor of $\pi/a$, also setting $k = 2p-1$, $l = 2q-1$, for $p,q \in \mathbb{N}$:
$$
F = \frac{16F_0}{\pi^2}\sum_{p=1}^\infty\sum_{q=1}^{\infty} X_{pq}(x)\sin\frac{(2p-1)\pi y}{a} \sin\frac{(2q-1)\pi z}{a}
$$
where 
$$X_{pq}(x) = - \frac{\cosh\frac{m\pi}{2} e^{m(x-1/2)\pi/a}}{(2p-1)(2q-1)\sinh(m\pi)}  +\frac{\cosh\frac{m\pi}{2}e^{-m(x-1/2)\pi/a}}{(2p-1)(2q-1)\sinh(m\pi)} $$
and
$$
 m =\sqrt{(2p-1)^2+(2q-1)^2}.
$$
