Your problem occurs when the upper limit of the integration exceeds k0. Are you sure the functions are defined correctly for that situation?
If the upper limit is taken as k0, the following code
%Enviroment Properties
f=15e9;
c=2.99792458e8;
lambda=c/f;
mu0=4*pi*1e-7;
eps0=1.0/c^2/mu0;
eta0=sqrt(mu0/eps0);
epsr=5;
omega=2*pi*f;
k0=omega*sqrt(mu0*eps0);
k2=omega*sqrt(mu0*eps0*epsr);
%square patch array properties
l=1.8e-3;
alpha_Ett=1.02*(l^3);
p=2e-3;
N=1/(p^2);
r=0.6956*p;
chieexx=(N*alpha_Ett)/(1-(N*alpha_Ett/(4*r)));
%source and observation points
z0=3e-3;
z=2e-3;
y=0;
%x = 5e-6:1e-4:5e-3;
x = logspace(-5.301,-2.301);
lo = 0; hi = k0; % Integration limits
for i = 1:numel(x)
pho=sqrt(x(i).^2+y.^2);
f1 = (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* (((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2))*(k2^2).*chieexx)+ ((k0^2)*(-2i+(sqrt(k0^2-kp.^2)).*chieexx)))))./((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2)).*(k2^2).*chieexx)+ ((k0^2).*(2i+(sqrt(k0^2-kp.^2)).*chieexx)))))) .* exp(1i.*(sqrt(k0^2-kp.^2)).*(z+z0)) );
f2= (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* exp(1i.*(sqrt(k0^2-kp.^2)).*(-z+z0)) );
I1= quadgk(f1,lo,hi);
I2= quadgk(f2,lo,hi);
E(i)= (eta0/(4*pi*k0)).*(abs(I1+I2));
hold on
end
semilogy(x,abs(E),'Or')
set(gca,'FontSize',14,'LineWidth',3)
xlabel('x [m]');ylabel('|Ez| [V/m]')
produces this result
