How to use integral2 when the integrand is a array?

(x,y) is a generated integrand as the following codes. d is a parameter in Tuu. I want to get a set of value of the integration "integral2(@Tuu,0,pi/2,0,pi/4)" with different values of d. And thus I set d=1e-9:1e-10:3e-9 as in the codes. However the codes "integral2(@Tuu,0,pi/2,0,pi/4)" gives the error "insufficient number of inputs". Why? How to solve this problem? Many thanks!

The codes of Tuu(x,y) are as following:

function U=Tuu(x,y)
d=1e-9:1e-10:3e-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

[Walter Roberson]

>> integral2(@Tuu,0,pi/2,0,pi/4)
Matrix dimensions must agree.
Error in Tuu (line 17)
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));

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

d=1e-9:1e-10:3e-9;
output = arrayfun(@(D) integral2(@(x,y) Tuu(x, y, D), 0,pi/2,0,pi/4,'reltol', 2e-4), d);

function U=Tuu(x, y, d)
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

If you try to use a smaller relative tolerance then you will get warning messages about using too many iterations. Your integrals are in the range of 2E20 so they do not converge well.

[Henan Fang]

@Walter Roberson Thank you very much!