Electric potential - Batygin, Toptygin book

I found in Batygin and Toptygin problems in electrodynamics a particular problem:

Density of charge is $$\rho = \rho_{0}\cos{\alpha x} \cos{\beta y} \cos{\gamma z}$$ in whole space. Find electric potential $$\phi$$.

Any hints where to start with? Direct integration isn't really pleasant thing to do by hand here. I can't think of any useful Gauss surface here neither.

Ok, I've found one way to get the right result but it requires a little cheating. First let's look at the Poisson's equation:

$$\Delta \phi = -4\pi\rho$$

My anzatz here is: potential $$\phi(x, y, z)$$ has to be a product: $$\phi(x, y, z) = \phi_{x}(x)\phi_{y}(y)\phi_{z}(z)$$ and that these potentials $$\phi_{i}(i) = A_{i}\cos{\xi_{i}x_{i}}$$,

where $$A_{i}$$ are constants, $$\xi_{i} = \{\alpha, \beta, \gamma \}$$ and $$x_{i} = {x, y, z}$$ for $$i = \{1, 2, 3\}.$$

I plugged it to Poisson's equation and got an answer:

$$\phi(x, y, z) = \frac{4\pi \rho_{0}}{\alpha^2 + \beta^2 + \gamma^2}\cos{\alpha x} \cos{\beta y} \cos{\gamma z}$$

Sadly, solution requires fortunate guess and I'm really curious if anybody knows more analytic approach.

You can solve the Poisson equation directly for this case using a Fourier transform. Using the convention $$\psi(\mathbf k) \equiv \mathfrak F[\phi] = \int_{\mathbb R^3}d^3 \mathbf x \ e^{-i \mathbf k.\mathbf x} \phi(\mathbf x)$$ for the 3D Fourier transform, you can take the Fourier transform of both sides of the Poisson equation to get: $$\mathfrak F[\nabla^2 \phi] = -4\pi \mathfrak F[\rho]$$ $$-|\mathbf k|^2 \psi(\mathbf k) = -4\pi \mathfrak F[\rho]$$ Which results in: $$\psi(\mathbf k) = \frac{4\pi}{|\mathbf k|^2} \mathfrak F[\rho] \quad \mathrm{(*)}$$ Now if we calculate $$\mathfrak F[\rho]$$, we can get $$\psi(\mathbf k)$$ from $$\mathrm{(*)}$$, and consequently find $$\phi(\mathbf x)$$ through an inverse Fourier transform. To calculate $$\mathfrak F[\rho]$$, note that for any function $$f: \mathbb R^3 \rightarrow \mathbb R$$ which can be written as a product of the form $$f(\mathbf x) \equiv f_x(x) f_y(y) f_z(z)$$, the 3D Fourier transform is simply the product of the 1D Fourier transform of each $$f_i$$ (why?).
Using this, you can easily show (exercise!) that the Fourier transform of $$\rho(\mathbf x) = \rho_0 \cos(\alpha x)\cos(\beta y)\cos(\gamma z)$$ is simply: $$\mathfrak F[\rho] = \pi^3\rho_0 \Big(\delta(k_x-\alpha)+\delta(k_x+\alpha)\Big)\Big(\delta(k_y-\beta)+\delta(k_y+\beta)\Big)\Big(\delta(k_z-\gamma)+\delta(k_z+\gamma)\Big)$$ Plugging this into $$\mathrm(*)$$ and taking the inverse Fourier transform gives: $$\phi(\mathbf x) = 4\pi^4 \rho_0 \int_{\mathbb R^3} \frac{d^3 \mathbf k}{(2\pi)^3} \ \frac{e^{i \mathbf k. \mathbf x}}{|\mathbf k|^2} \Big(\delta(k_x-\alpha)+\delta(k_x+\alpha)\Big)\Big(\delta(k_y-\beta)+\delta(k_y+\beta)\Big)\Big(\delta(k_z-\gamma)+\delta(k_z+\gamma)\Big)$$ Calculating this integral using the properties of the delta distribution, the final result is: $$\phi(\mathbf x) = \frac{4 \pi \rho_0}{\alpha^2+\beta^2+\gamma^2}\cos(\alpha x)\cos(\beta y)\cos(\gamma z)$$ as derived through your own method. I've left the final calculations for you to go through yourself.