# Numerical solution of the BCS gap equation, with Coulomb potential

I'm interested in the excitons condensation, which is described by an equation which is very similar to the standard BCS gap equation. I will be referring to this article: http://arxiv.org/abs/cond-mat/9501011. The only differences w.r.t. the standard BCS theory are:

1. The momentum dependency of the gap cannot be neglected, so we write $\Delta = \Delta(k)$.
2. The potential must be model realistically, in this case after a Coulomb potential.

In the end the gap equation reads:

$$\Delta (k) = \int d^2 k' \underbrace{\frac{e^{-d \left| k - k' \right|}}{\left| k - k' \right|}}_{V_{k k'}} \frac{\Delta(k')}{\sqrt{\Delta^2 + (k'^2 - \mu)^2}}$$

Knowing all the the parameters ($d$, $\mu$), clearly this equations needs to be solved numerically for $\Delta (k)$.

My first guess (and this approach does work, indeed, in many similar cases) would be of choosing a "trial" $\Delta$, iterating over and over, hoping sooner or later to converge to a solution which is not the banal $\Delta=0$.

This approach does not work here, because the potential $V_{k k'}$ has singularities. The article I am referring to is quite cryptical about it. They mention an unpublished paper as far as the numerical calculations are involved, and simply say that "the equation is converted in matrix form", as:

$$\Delta_i = M_{i j} \Delta_j$$

as expected, and then "the logarithmic singularities are treated separately". Clearly the singularities are on the diagonal of the matrix when the equation is discretized. This process is completely obscure to me. I suppose I cannot just separe the divergent diagonal part and the finite part. I tried parametrizing the divergency as $\frac{1}{\epsilon}$ and letting $\epsilon$ go to zero at the end of the calculation but to no avail. Is there some standard techniques to numerically solve my equation?