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.