Also shown in Fig.3,4 are the numerical results for the full
self-consistent integral equation, (14). This we have solved
iteratively after the fashion of Broyles [22]: an initial
guess input function 
is used in computing the right hand side of Eq.(14) giving an
output function 
, and the next input function is chosen to be:
where 
is a constant chosen between 0 and 1 for
optimal convergence of the process. For computational purposes,
the function 
has to be specified on a two
dimensional grid 
(assuming unbroken symmetry about
the axis normal to the plates), with h being a fixed parameter
during the iterations. The function is obtained by fitting bi-cubic
splines between grid points.
A fine grid means more grid points where evaluations of the function
n has to be made, whereas a coarse one would not supply a
satisfactory accuracy (which we define to be convergence of the
force to within 
of the maximum value
of 
for any given force curve at fixed 
).
The bottleneck of the iterative
process is the integral contained in 
,
which is of high dimension, and makes precise evaluations very difficult.
This problem is alleviated by splitting the integral
into two more straightforward ones:
where 
, 
,
and f and g are two dimensionless
functions:
Substituting f and g into Eq.(34), together with the relations
between 
's, gives the original Eq.(15).
Note that f(t,t') can be calculated and tabulated beforehand to speed up the iteration
procedure, and 
is a simple integral
with its integrand specifying a straight line segment in the
plane.
It turns out the number of iterations needed to
converge Eq.(33) is of the order 5-15 with
between 0.1 and 0.5, and rather
insensitive to the initial guess function , which is therefore
always set to the constant bulk value 
. Not surprisingly,
the process is a lot faster for smaller concentrations ( 
),
and it gets increasingly difficult when approaches 2, which we
found to be the maximum value for which converged results could be obtained
in reasonable time (5-10 hrs CPU per value on a DEC
Alpha 600/5/266). In order to have
the fastest convergence, a two dimensional adaptive grid was adopted
which becomes finer
wherever the function n varies most rapidly.
After numerically solving Eq.(14) for ,
and following Eq.(4), we can find:
We found that this leads, with a high degree of numerical accuracy and
irrespective of the value of to:
so that the resulting pressure 
is once more identical
to the second virial result [21, 18], which is known to be exact in the
limit of 
which we have adopted.
The depletion force for general separation h is then given by equations (3, 37) as:
As in the perturbative treatments of Section 3, the force between spheres of radius R is then found by the Derjaguin integral, Eq.(1); the interaction energy between two such spheres can then be found by a further integration.
