If I understand what you're trying to do, you should note that FourierTransform[f[x,y,z],{x,y,z},{kx,ky,kz}] does return a function of kx,ky, and kz, even if it remains in that form. That is, running the following code
Ff[kx_, ky_, kz_] :=
FourierTransform[f[x, y, z], {x, y, z}, {kx, ky, kz}];
InverseFourierTransform[Ff[kx, ky, kz], {kx, ky, kz}, {x, y, z}]
correctly returns f[x,y,z]. I'd presume then you can work with the Fourier transform symbolically, but it seems more practical to just use DSolve or NDSolve if you're solving a system of PDEs.
