2d Chebyshev quadrature for approximation of unsolvable integral

Hey :) I am trying to approximate a double integral via Chebyshev-Gauss. My code is:

if 
 hbar = 6.582119e-01; % Planck constant (meV ps)
Delta = 1.0; tau = 0.5; o = 3./hbar; A = 1; p = 2.*sqrt(log(2))./tau;
n = 5; m = 5;
[omega,d] = meshgrid(2:0.25:10,0:0.25:10);
syms k l
x_k = 1./2.*cos((2.*k-1)./(2.*n).*pi)+1./2;
x_l = 1./2.*cos((2.*l-1)./(2.*m).*pi)+1./2;
w_1 = pi./n;
w_2 = pi./m;
f = (exp(-((o-omega-x_k-x_l)./(2.*p)).^2).*sin(d.*(omega+x_k+x_l))-exp(-((o-omega+x_k+x_l)./(2.*p)).^2).*sin(d.*(omega-x_k-x_l))).*((1-x_k.*x_l).*Delta.^2)./(x_k.^2.*x_l.^2);
F = symsum(f,k,1,m);
G = tau./(8.*omega).*A.*sqrt(pi./log(2)).*1./4.*w_1.*w_2.*symsum(G,l,1,n);
plot3(d,omega,G);
end

*the actual function i want to integrate is

if 
1./((1-x.^2).*(1-y.^2)).*f
end

though. I set the Chebyshev variables (x,y) to 0.5x+0.5 and 0.5y+0.5 because I want to integrate from 0 to 1 (and the Chebyshev approximation is defined for integrations from -1 to 1). The goal is to approximate the integral and then plot my function depending on d and omega... I would appreciate any kind of help...(or maybe also different approximation attempts if you have...) I am just desperate with the calculation of this integral. Edit: maybe I should have said: I don't get any error, matlab is just compiling and compiling but nothing happens and after some time my computer just shuts down...