A macroscopic model to describe the dynamics of ion transport in ion channels
is the Poisson-Nernst-Planck(PNP) equations. In this paper, we develop a
finite-difference method for solving PNP equations, which is second-order
accurate in both space and time. We use the physical parameters specifically
suited toward the modelling of ion channels. We present a simple iterative
scheme to solve the system of nonlinear equations resulting from discretizing
the equations implicitly in time, which is demonstrated to converge in a few
iterations. We place emphasis on ensuring numerical methods to have the same
physical properties that the PNP equations themselves also possess, namely
conservation of total ions and correct rates of energy dissipation. We describe
in detail an approach to derive a finite-difference method that preserves the
total concentration of ions exactly in time. Further, we illustrate that, using
realistic values of the physical parameters, the conservation property is
critical in obtaining correct numerical solutions over long time scales
