Warning: The maximum number of function executions has been reached (10000). The results failed the global error test.

If you approach the singularity(z=0) during double integration, the following error occurs.

Warning: The maximum number of function executions has been reached (10000). The results failed the global error test.(It is a translation of the expression)

a picture is attaced to help understanding.

I wonder why the error occurs and how to fix it.

Below is the code I used.

Tau1=1;  
partition_z=100000;  
rmax=1;  
zmax=1; 
qbz=zeros(1,partition_z+1);
thetamin=0;
thetamax=pi;
Rmin=0;
Rmax=1;
r=0;
parfor k=1:1:partition_z+1
    z=zmax*(k-1)/partition_z;
    funb_z_withoutconstant=@(theta,R) (2/pi)*(exp(-sqrt(r.^2-2*r.*R.*cos(theta)+R.^2+z.^2).*Tau1).*R.*z.^2)./((sqrt(r.^2-2*r.*R.*cos(theta)+R.^2+z.^2)).^4);
    qbz(k)=integral2(funb_z_withoutconstant,thetamin,thetamax,Rmin,Rmax);
end

z=0:zmax/partition_z:zmax;
plot(z,qbz)
hold on
title( 'nomalized qbz')
xlabel('z')