$\newcommand{\P}{\mathbb{P}}$ I have a $2$-D stationary Fokker-Planck equation $$\frac{\partial^2 \P(A,B)}{\partial A^2}+\frac{\partial^2 \P(A,B)}{\partial B^2}=\frac{\partial f_1(A,B) \P(A,B)}{\partial A}+\frac{\partial f_2(A,B) \P(A,B)}{\partial B},$$ where $f_1$ and $f_2$ are polynomials.
I would like to compute the solution with the normalization condition $\int \P(A,B)\mathrm{d}A\mathrm{d}B=1$. I have two question:
Is this possible without any boundary conditions?
If I chose to use finite differences method, a very naive choice would be to do something like this:
$$\frac{\P_{(a+1,b)}-2\P_{(a,b)}+\P_{(a-1,b)}}{\Delta A^2}+\frac{\P_{(a,b+1)}-2\P_{(a,b)}+\P_{(a,b-1)}}{\Delta B^2}\\=\frac{\P_{(a+1,b)}f_{1(a+1,b)}-\P_{(a-1,b)}f_{1(a-1,b)}}{2\Delta A}+\frac{\P_{(a,b+1)}f_{2(a,b+1)}-\P_{(a,b-1)}f_{2(a,b-1)}}{2\Delta B},$$
where $(a,b)$ are indexes of 2D uniform mesh, with an extra equation given by the discretization of the normalization condition.
Is this a good way to solve the problem ?
Your equation may be rewritten in somewhat more standard terms by viewing it as a problem on $\mathbb{R}^2$. For this, just compute the derivatives on the left-hand side: $$ \Delta p = \partial_A (f_1 p) + \partial_B (f_2 p) \\ = (\partial_A f_1) p + f_1 \partial_A p + (\partial_B f_2) p + f_2 \partial_B p \\ = \left[f_1 \partial_A p + f_2 \partial_B p\right] + \left[(\partial_A f_1 + \partial_B f_2) p\right]. $$ Introducing $F = (f_1, f_2)^T$, you may rewrite the equation as $$ -\Delta p + F \cdot (\nabla p) + (\nabla \cdot F) p = 0 $$
This is a stationary convection-diffusion-reaction equation. In theory, you don't have boundary conditions for $p$, however the normalisation condition implies that $p$ will decay towards $0$ for $A, B \rightarrow \pm \infty$. So an approach that could work is to compute the solution on a "sufficiently large" box $$B = (-L, L) \times (-L, L) \subset \mathbb{R}^2, \quad L > 0,$$ where you impose homogeneous Dirichlet boundary conditions, i.e. $p|_{\partial B} = 0$.
Proving existence and uniqueness in the above setting depends on the polynomials $f_1$ and $f_2$, but I assume in general, the solution will only be unique modulo scaling by a constant. Although the normalisation constraint could be incorporated in a Finite Element framework (using Mixed Finite Elements more specifically), I think that it might be worth a try to simply apply the Finite Difference Method to the equation with Dirichlet boundary conditions and afterwards compute the integral $$C = \int\limits_B \! \mathrm{d}x~p$$ numerically and then scale your obtained solution by $C^{-1}$.
For specific choices of polynomials $f_1$ and $f_2$ there might even be simpler solutions to this problem, e.g. in the case that $\nabla \cdot F = 0$.