A hybrid of two numerical schemes is used.
The first scheme relies on the observation that because the recursion relations 9 and 10 are linear, the values of the flux and potential at any facet is a linear function of the the values at the first facet:
Here each of the 
, 
,
and 
are complicated functions of 
which can be numerically computed
for a particular configuration of facets using 9
and 10 by the following stratagem:
Let 
and 
be unknowns, but calculate from the
recursion
relations 9 and 10, 
and 
for the case that 
and 
;
and 
.
The 
can be calculated by setting 
and
and the 
are calculated similarly.
Because the final endpoint of segment 
is the initial point of
segment 
, 
and
(equation 8, the unknowns 
and
can be determined directly from the solution of the two linear
equations 24.
All 
, 
, and 
can then be determined
by recursion equations 9- 13.
While this method lends itself to rapid computation, it is numerically unstable when there is even a mildly disparate range of facet lengths. Truncation errors currently render it useless for more than around 20 facets. The problem is alleviated by using structures which can manipulate numbers with arbitrary precision, but the computations rapidly become slower as the number of facets in a surface increases and the necessary precision increases. Because of the lack of speed, this method is only used when the following scheme becomes unreliable.
The second method of solution is a conjugate gradient scheme[11]
which is typically an 
iterative scheme.
To implement the conjugate gradient, the recursion equations can
be written as a linear system:
where,
Conjugate gradient iteration works by squaring the difference of equation 25 and iterating the positive definite system by choosing an optimal set of iteration directions. Squaring the system also squares the condition number of the matrix, which makes the numerical precision worse for small condition numbers. Experience has shown that the stability of the conjugate gradient methods is increasingly worse as the length scales of a particular surface become more disparate, and when this occurs the first scheme is applied using arbitrary precision numbers.
Integration in time is performed explicitly with time steps being chosen dynamically so that no pair of edges may overrun an intervening edge and so that no edge velocity may change too rapidly.