  Electrostatic correlations and fluctuations in ionic systems can be described
within an extended Poisson-Boltzmann theory using a Gaussian variational form.
The resulting equations are challenging to solve because they require the
solution of a non-linear partial differential equation for the pair correlation
function. This has limited existing studies to simple approximations or to
one-dimensional geometries. In this paper we show that the numerical solution
of the equations is greatly simplified by the use of selective inversion of a
finite difference operator which occurs in the theory. This selective inversion
preserves the sparse structure of the problem and leads to substantial savings
in computer effort. In one and two dimensions further simplifications are made
by using a mixture of selective inversion and Fourier techniques.
