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}. $$