  In this work we design and analyze a free energy satisfying finite difference
method for solving Poisson-Nernst-Planck equations in a bounded domain. The
algorithm is of second order in space, with numerical solutions satisfying all
three desired properties: i) mass conservation, ii) positivity preserving, and
iii) free energy satisfying in the sense that these schemes satisfy a discrete
free energy dissipation inequality. These ensure that the computed solution is
a probability density, and the schemes are energy stable and preserve the
equilibrium solutions. Both one and two-dimensional numerical results are
provided to demonstrate the good qualities of the algorithm, as well as effects
of relative size of the data given.
