Help solving an integral numerically

I am trying to solve an integral where there is only an integration variable (kx in the code below) but it has to be solved in a 2D space defined by x and y coordinates. It corresponds to a Green's function in 2D.

The following is my approach defining it as a handle function, but it seems to take for ever to run.

clc; clear all;

f = 8500;
c0 = 343;
rho = 1.225;
omega = 2*pi*f;
k = omega/c0;
Z = -426;

lx = 0.1;
ly = 0.1;
x0 = 0;
y0 = ly/1000;
nx = 50;
ny = nx/2;
x = linspace(-lx,lx,nx);
y = linspace(0,ly,ny);

G = zeros(length(x),length(y));

integrand = @(x,y,kx) ((exp(1i*(kx*x + sqrt(k.^2 - kx.^2).*abs(y))))./sqrt(k.^2 - kx.^2));
Gz = @(x,y) integral(@(kx)integrand(x,y,kx), 0, 200, 'ArrayValued', true);

for ii = 1:length(x)
    for jj = 1:length(y)
            
        G(ii,jj) = (1j/4/pi)*Gz(x(ii),y(jj));
       
    end
end

Any help will be appreciated.

[Torsten]

Remove the 'ArrayValued',true - option.

[smonos]

If I remove it I get the following error:

Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued'

option to true.

[Torsten]

Then use

integrand1 = @(x,y,kx) real(((exp(1i*(kx*x + sqrt(k.^2 - kx.^2).*abs(y))))./sqrt(k.^2 - kx.^2)));
integrand2 = @(x,y,kx) imag(((exp(1i*(kx*x + sqrt(k.^2 - kx.^2).*abs(y))))./sqrt(k.^2 - kx.^2)));
Gz = @(x,y) integral(@(kx)integrand1(x,y,kx), 0, 200) + 1i*...
            integral(@(kx)integrand2(x,y,kx), 0, 200);

[smonos]

Many thanks!