Numerical integration of divergent functions

Question

I am having trouble with the numerical integration of a divergent function. For example, \begin{equation} n= \int f(x)\,dx = \displaystyle\int \dfrac{\Theta(x-\varepsilon)\,dx}{\sqrt{x-\varepsilon}} \end{equation} is not behaving good when calculated numerically: \begin{equation} n_\text{numerical} \approx f(x_j)(x_j-\varepsilon) + \sum_{k=j+1}^N f(x_k)\Delta x \end{equation} with \begin{equation} x_i = x_\text{min} + i\Delta x \end{equation} and $x_j$ being the smallest point larger than $\varepsilon$. Due to the divergent behavior of $f(x)$ near $\varepsilon$, the numerical result has discontinuities and jumps at every $\varepsilon=x_i$ when I sweep $\varepsilon$. The smoothness of this integration is crucial in the convergence of my self-consistent simulations. Can anybody suggest a simple trick to solve the problem?