Try :
integral2(@Tuu,0,pi/2,0,pi/4) % calling if the function
function U=Tuu(x,y)
d=2e-9;
mu=8;
delta=10;
vh=16;
HBAR=1.05457266e-34;
ME=9.1093897e-31;
ELEC=1.60217733e-19;
Kh=2.116e10;
kc=sqrt(2.*ME.*ELEC./HBAR.^2);
ku=kc.*sqrt(mu+delta);
kd=kc.*sqrt(mu-delta);
puu1=sqrt(ku.^2-ku.^2.*sin(x).^2+kc.^2.*vh);
puu2=sqrt(ku.^2-ku.^2.*sin(x).^2-kc.^2.*vh);
quu1=sqrt(ku.^2-(Kh-ku.*sin(x).*cos(y)).^2-ku.^2.*sin(x).^2.*sin(y).^2+kc.^2.*vh);
quu2=sqrt(ku.^2-(Kh-ku.*sin(x).*cos(y)).^2-ku.^2.*sin(x).^2.*sin(y).^2-kc.^2.*vh);
U=1/4.*ku.*sin(x).*(real(puu1).*exp(-2.*imag(puu1).*d)+real(puu2).*exp(-2.*imag(puu2).*d)+real(quu1).*exp(-2.*imag(quu1).*d)+real(quu2).*exp(-2.*imag(quu2).*d)+((real(puu1)+real(puu2)).*cos((real(puu1)-real(puu2)).*d)-(imag(puu1)-imag(puu2)).*sin((real(puu1)-real(puu2)).*d)).*exp(-(imag(puu1)+imag(puu2)).*d)-((real(quu1)+real(quu2)).*cos((real(quu1)-real(quu2)).*d)-(real(quu1)-imag(quu2)).*sin((real(quu1)-real(quu2)).*d)).*exp(-(imag(quu1)+imag(quu2)).*d));
end
