Solve for F=(F1;F2;F3)

Where F1,x + F2,y +F3,z =f(x,y,z) = gaussian

I doubt the answer is unique. You must have some symmetry conditions.

It is a 3D Gaussian.

Maybe express everything in spherical coordinates, then it will simplify.

Of course, the system is underdetermined. Suppose you have already a solution of div(F) = f, then add the curl of an arbitrary vector field A, i.e., set F'=F+curl A. So you have div(F') = div(F)+div curl A = div(F) = f. This will not be eliminated by integrability conditions. Boundary conditions may render the solution unique.

According to Helmholtz decomposition in 3 dimensions, you can allways write F= -grad (u)+curl (A), where u is a scalar field and A is a vector field. So, operating with the divergence in both sides you get lap (u)=-f, where "lap" stands for Laplacian. Then you use boundary conditions on F to select those for u and A. For instance, in the relevant case of full 3-dimensional space demanding that F and all its derivatives vanish at infinity, one solution to your problem is simply F=-grad (u), where u is the convolution of your Gaussian with 1/|r|. Additional solutions can be found by adding to the above solution the curl of any vector field A so that curl (A) vanishes at infinity (i'm not sure whether such a vector field A exists, if we additionally demand it to be smooth).

Thanks for these comments. Of course the solution is not unique and boundary conditions are needed. I want the solution F and its derivatives to vanishe at infinity.

I have now a partial answer to this question.

Let us denote Phi(x)=int_{-infty}^x exp(-t^2) dt.

If we make a change of variables so that

f(x)=exp(−sumi xi2)

then taking Fi(x) = 1/n Phi(x_i) exp( - sum_{ji} xj2)

gives div(F) = f.