I want to calculate numerically the triple integral of a function which is defined itself by the product of two integrals. I managed to get a result with Matlab after a few minutes but I did not manage with Maple even after several hours of calculation. Here is the script.

pmint := proc(x1, y1, z1) local nu_0, nu_m, Qm, pm; nu_0 := sqrt(k^2 - k_0^2); nu_m := sqrt(k^2 - k_m^2); Qm := omega*rho_me*I; pm := 2*Qm*rho_0*(1 + tanh(nu_m*d))*exp(-nu_0*zs)*exp(-nu_m*d)*cosh(nu_m*(z + d))*BesselJ(0, k*sqrt((xs - x)^2 + (ys - y)^2))*k/(rho_me*nu_0 + rho_0*nu_m*tanh(nu_m*d)); `~`[evalf@Int]((Re + Im*I)(eval(pm, [k_0 = 0.018371886863098, xs = 0, ys = 0, zs = 2, rho_0 = 1.2, rho_me = 1.528516816439260 + (-1)*1235.297048886680*I, k_m = 0.490806242885258 + (-1)*0.490314205803914*I, d = 0.05, omega = 2*Pi, x = x1, y = y1, z = z1])), k = 0 .. infinity, epsilon = 0.3*10^(-5), method = _d01amc); end proc;
dom1 := [x1 = -2/2 .. 2/2, y1 = -2/2 .. 2/2, z1 = -0.05 .. 0];
evalf(Int(pmint(x1, y1, z1)*conjugate(pmint(x1, y1, z1)), dom1, digits = 6));

I get the following value with Matlab : 9.8202895

with a relative Tolerance of 1e-6 and an absolute tolerance of 1e-12 in the integral and integral3 Matlab routines..

Could you help me ?

Thanks