  Solving the stationary nonlinear Fokker-Planck equations is important in
applications and examples include the Poisson-Boltzmann equation and the two
layer neural networks. Making use of the connection between the interacting
particle systems and the nonlinear Fokker-Planck equations, we propose to solve
the stationary solution by sampling from the $N$-body Gibbs distribution. This
avoids simulation of the $N$-body system for long time and more importantly
such a method can avoid the requirement of uniform propagation of chaos from
direct simulation of the particle systems. We establish the convergence of the
Gibbs measure to the stationary solution when the interaction kernel is bounded
(not necessarily continuous) and the temperature is not very small. Numerical
experiments are performed for the Poisson-Boltzmann equations and the two-layer
neural networks to validate the method and the theory.
