" Arrays have incompatible sizes for this operation" when computing nested integral

function output = tsp_up(Z)     
ap = active_p(Z); % A constant
disp(ap);
theta = const.alpha / 2;
fun = @(x,y) 1 ./ (1 + const.gamma_B ^(-1) .* (x.^(theta)).* (y.^(- theta))) .* exp(- ((Z - const.gamma_B).* const.Noise .* (x.^(theta))) ./ (const.gamma_B .* const.P));
tsp_r = @(y) y .* exp( -pi * ap * const.lambda_BS .* integral(@(x)fun(x,y),0,Inf) - pi .* y.^2 .* const.lambda_BS);
result = integral(tsp_r,0.01,100);
output = 2 * pi * const.lambda_BS * result;
end

y = tsp_up(5.1);

Here I want to compute the nested Integral with variable x(from 0 to Inf) and y(from 0.01 to 100). Everything else is constant. The " Arrays have incompatible sizes for this operation" error has occured in line 5(fun) . I can't see what it has anything to do with the array size. Does anyone know where I am wrong? Any kind of advice is sincerely appreciated!

I try to modify the domain of x and y to be (0.01 to 100) both, this time the error disappear but I totally don't understand what is happening.

[Walter Roberson]

I recently explained this at https://www.mathworks.com/matlabcentral/answers/882478-integral-inside-integral-is-taking-too-much-time#comment_1687617

[chenjie ruan]

Thank you so much! It works after modify it as

result = integral(tsp_r,0.001,100,'ArrayValued', false);

But I am still confused about why it won't work when the outer and inter integral is both arrayvalued. Is it because the inner integral does not have a close-form so we have to let it compute the value at every y? I really appreciate it if you could answer this.

[Walter Roberson]

Unfortunately you did not happen to include enough for us to be able to test the code.

When you use integral() it always does numeric integration, so closed form or not does not matter.

[chenjie ruan]

Huge thanks for your answer!