Numerical double integral with singularities

Dear MATLAB users,

This is a mathematical question which I could not solve it properly. So I come for your help and any suggestions would be welcomed.

My problem has be simplified to this point that I have to finish the following double integral using numerical methods for given parameter ξ (and ).

I tried three functions to finish this integral but I fail to make it. I received many warnings.

Here are my codes. What I want to do is plot the function vs ξ. The problem is that I don't know how to deal with the singularities of the integrand.

xi = -3:0.1:3;
n = numel(xi);
fvals = zeros(1, n);
for i = 1:n
    xi0 = xi(i);
    fun = @(x, y) 1./(sqrt(1-x.^2).*sqrt(1-y.^2).*sqrt(1-(xi0-x-y).^2));
    fvals(i) = quad2d(fun, -1, +1, -1, +1);
    %fvals(i) = integral2(fun, -1, +1, -1, +1);
    %fvals(i) = dblquad(fun,-1,+1,-1,+1);
end
figure; plot(xi, fvals, '.-')

Here are the warnings.

...

> In quad2d (line 244)

警告: 已达到函数计算的最大数目(2000)。结果将使全局误差测试失败。

> In quad2d (line 244)

警告: 已达到函数计算的最大数目(2000)。结果将使全局误差测试失败。

> In quad2d (line 244)

警告: 已达到函数计算的最大数目(2000)。结果将使全局误差测试失败。

> In quad2d (line 244)

警告: 已达到函数计算的最大数目(2000)。结果将使全局误差测试失败。

> In quad2d (line 244)

警告: 复数 X 和/或 Y 参数的虚部已忽略

...

And here are the resutls.

Best regards

Qilin.