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