3D integral failure

I am trying to sovle the following integration but integral3 keeps failing. Is there any work around for this. The expression for integration is:

with

where Tr[...] represent trace of ... matrices multiplication, b = 1/(kB . Q), kB, hbar, Q are constants, and

and

with h(x,y) is a 3-by-3 matrix given in below code F.m. I call integral3() as shown in runner.m code.

%runner.m
clear; clc;

% parameters
Q = 34.7;
eta = 0.025;
hbar = 6.582119569e-13;

% limits of int
xmin = -2*pi/3; xmax = 4*pi/3;
ymin = -2*pi/sqrt(3); ymax = 2*pi/sqrt(3);
zmin = 3; zmax = 15;

fun = @(x,y,z) hbar/((2*pi)^3*Q) * F(x,y,z,Q,eta);
I = integral3(@(x,y,z)arrayfun(fun,x,y,z),xmin,xmax,ymin,ymax,zmin,zmax);

It gives that the result is NaN.

%F.m
function F = F(x,y,z,Q,eta)
%=============== constants ===========
J = 1;
S = 1;
Delta = 2;
DbyJ = 2/sqrt(3);
hbar = 6.582119569e-13;
kB = 8.617333262145e-2;
b = 1/(kB*Q);


% ================ h(x,y),vx,vy ============
%some parameter for h(x,y)
t1 = sqrt(3)*x;
t2 = (3*y+t1)/2;
D = 1+1i*DbyJ;
fk1 = 1+exp(1i*t1);
fk2 = 1+exp(1i*t2);
fk1k2 = 1+exp(1i*(t1-t2));

%h(x,y)
h = J*S*[4*Delta, -D*conj(fk1), -conj(D*fk2);
    -conj(D)*fk1,4*Delta, -D*fk1k2;
    -D*fk2, -conj(D*fk1k2),4*Delta];

%vx, vy calculated using syms expression
vx = 1/hbar*[                                                                    0,                                   3^(1/2)*J*S*exp(-3^(1/2)*conj(x)*1i)*(1 + DbyJ*1i)*1i, -(3^(1/2)*J*S*exp(- (conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*1i)/2
    3^(1/2)*J*S*exp(3^(1/2)*x*1i)*(- 1 + conj(DbyJ)*1i)*1i,                                                                                        0,                     -(3^(1/2)*J*S*exp(- (y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*1i)/2
    -(3^(1/2)*J*S*exp((y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*1i)/2, -(3^(1/2)*J*S*exp((conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*1i)/2,                                                                                          0];

vy = 1/hbar*[                                                            0,                                                                               0, -(J*S*exp(- (conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*3i)/2
    0,                                                                               0,                      (J*S*exp(- (y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*3i)/2
    -(J*S*exp((y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*3i)/2, (J*S*exp((conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*3i)/2,                                                                                  0];

% ===================== GA, GR ================
GR = inv((z+1i*eta)*eye(3)-h);
GA = GR';
G = GA-GR;
term1 = b*z/(4*(sinh(b*z/2))^2);
F = -trace((z*vy)*G*vx*G)*term1;