How to numerically integrate bessely function

Hi All, i want to get the integral result for bessely function as

int{a}{b}int{a}{b}bessely(1, x1-x2)dx1dx2

but when numerically calculate that using Gaussian quadrature points, the values in bessely(1, 0) is a infinite number, and the limit as limit{x->0}bessely(1, x) are not existed. But when i used the quad2d to calculate this integral, a correct answer can be obtained. So i like to know how the results can be numerically obtained? Thanks.

My code is as

clear
clc
close all
a=-8;
b = 6;
% covpara=1;

%%number of quadrature point
n=50;

% % % % the whole gaussian quadrature point
[Q, W]=LegGaus(n,a,b);

covfun_A=@(x1,x2) abs(x1-x2).*bessely(1,abs(x1-x2));

Matrix_A = bsxfun(covfun_A,Q.',Q); % f = @(x,y)....

%A_1 impose the diagnoal elements as a very small number, i.e., bessely 
Matrix_A(1:length(Matrix_A)+1:length(Matrix_A)*length(Matrix_A))=0; 

%A_1 is the numerical result using Gaussian quadrature, using tensor
%product of 1D to get 2D result
A_1 =(W)'*Matrix_A*(W);

%A_2 is the 'exact' result from matlab
A_2= (quad2d(@(x1,x2) (abs(x1-x2).*bessely(1,abs(x1-x2))),a,b,@(x1)x1,b) +...
    quad2d(@(x1,x2) (abs(x1-x2).*bessely(1,abs(x1-x2))),a,b,a,@(x1)x1));

A_2/A_1

LegGaus is a sub function in getting the weight and quadrature points. If you need that, i can mail it to you.

thanks.

[David Zhang]

BTW, i dont know how to choose Matrix_A(1:length(Matrix_A)+1:length(Matrix_A)*length(Matrix_A))=0;

if i increase the number of Gaussian points to a very large number, these two results can be quite the same, however, it is time-consuming.

[Walter Roberson]

Context: http://www.mathworks.com/matlabcentral/answers/27994-why-dblquad-can-not-be-used-to-evaluate-the-bessel-function-of-the-second-kind-bessely