Recursive algorithm for numeric integration with singularities

Question

As part of my research, I have to integrate Green's Function to get charge density but Green's function has integrable-singularity near band-edge and integration has been a headache with normal trapezoidal/simpson rule. I have never done any formal course or informal training on numerical techniques except error analysis like butcher tableau.

Based on the internet search, typical example of $\int_0^1\frac{1}{\sqrt{\sin(x)}}dx$ is given where we know the point of discontinuity and then by transformation of variable this can be tackled.

More generally, we don't know the location of integrable singularity (like in Green's function) and then it becomes difficult to use such a transformation.

So far what I have been doing is to increase the number of sampling points but this slows down the code significantly. What would really help is if the sampling increases near the point of discontinuity so I wrote a recursive code and tried it on the $\frac{1}{\sqrt{x-0.5005}}$ example.

The recursion works as follows:

The function that perform integration is called with end limits, first trapezoidal integration is performed with just two end points and if the integral is greater than tolerance, two more calls are made with the limits change to cover first and second half of the x-axis that was supplied to the original function, this is expressed in the lines:

function [I,x] = integrate(x1, x2, tol)
    persistent n
    if isempty(n)
        n = 0;
    end
    n = n+1;

dx = x2-x1;
I = (f(x1)+f(x2))*dx/2;
if(I<tol)
else
    I = integrate(x1, x1+dx/2, tol)+integrate(x2-dx/2, x2, tol);
end
x = n;
end

I compare the result from recursion with the simpson (from matlab exchange) and trapezoidal method and the entire Matlab script is as follows:

tol = 1e-6;
[I,n] = integrate(0,1, tol);

format long;
actual = 2*sqrt(0.4995);
fprintf("Exact integral %.10f vs recursive integral %.10f \n", actual, I)
x = linspace(0,1,n+1);
f_sim = f(x);
I_sim = simps(x', f_sim', 1);
f_trapz = cumtrapz(f_sim/n);
I_trapz = f_trapz(end);

fprintf("Exact integral %.10f vs simpson integral %.10f \n", actual, I_sim)
fprintf("Exact integral %.10f vs trapezoid integral %.10f \n", actual, I_trapz)

e1 = abs(I-actual)/actual;
e2 = abs(I_sim - actual)/actual;
e3 = abs(I_trapz - actual)/actual;
fprintf("Recursion(%.03e) vs Simpson(%.03e) vs Trapz(%.03e) error  \n", e1, e2, e3)

function [I,x] = integrate(x1, x2, tol)
    persistent n
    if isempty(n)
        n = 0;
    end
    n = n+1;

dx = x2-x1;
I = (f(x1)+f(x2))*dx/2;
if(I<tol)
else
    I = integrate(x1, x1+dx/2, tol)+integrate(x2-dx/2, x2, tol);
end
x = n;
end

function fun = f(x)
fun = 1./sqrt(x-0.5005).*(x>=0.5005);
end

The output that i get is as follows:

Exact integral 1.4135062787 vs recursive integral 1.4135064792 
Exact integral 1.4135062787 vs simpson integral 1.4137763950 
Exact integral 1.4135062787 vs trapezoid integral 1.4135491522 
Recursion(1.418e-07) vs Simpson(1.911e-04) vs Trapz(3.033e-05) error  

So for same number of function calls, recursion is indeed giving far greater accuracy than other two methods. All methods use 4 million points.

Can someone point to resources that can help me write or get or understand such recursive codes for integration with more optimization and for huge matrices?