Double integration of a function formed by multiplying two function handles

The code involves first the integration of two functions with respect to z and r1 respectively as shown below.

func_L_Iagg = @(x,y,z) (-78.709.*((z.^-0.5 + z.^-0.333)/(1+0.3162.*z.*(x.^-1+y.^-1))));
L_Iagg_x_y = @(x,y) exp(integral(func_L_Iagg,y,inf));


func_L_Ir1 = @(x,y,r1) ((0.667)./((r1.^0.333).*(x.^0.667).*(1+((3.1623)./((x.^-1+y.^-1).*r1)))));
L_Ir1_x_y = @(x,y) integral(func_L_Ir1,0,x);

As per my understanding of the meaning of function handle, the statements above will generate two function handles. Next part of the problem is to create a function by multiplying the above two generated results and double integrate it over the given limits.

xmin = 0;
xmax = inf;
ymin = @(x) x;  
ymax = inf;
fun = @(x,y) L_Ir1_x_y(x,y)*L_Iagg_x_y(x,y);
C_BK_FS = integral2(fun,xmin,xmax,ymin,ymax);

Error: A and B must be floating-point scalars.

May be my understanding of the term function handle is wrong or I am using the integral command in a wrong way. Can someone please help me out with this?

Thanks in advance!

Answer 1

[Walter Roberson]

The call to integral2 itself is not wrong. However, integral2(),

The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations.

You take the x and y data and use it as bounds in integral(), but it is non-scalar values and you cannot use non-scalar values as upper or lower bounds for integral() or integral2()

Your two individual functions need to arrayfun() over the bounds. And your multiplication of the two functions will need to be .* instead of *