I posted an answer to your deleted question.
I think that the moderator who deleted it was wrong because:
1. It was answered.
2. The content was distinct from the previous question. Not to mention that the old thread is difficult to follow (and loads slowly) being very long.

So, here is the answer again.

f:=exp(-(sqrt(4*x^2+4*y^2+4*z^2)^3));
F:=eval(f, [x=r*sin(v)*cos(u), y=r*sin(v)*sin(u), z=r*cos(v)]);
int(F*r^2*sin(v), r=0..2, u=0..2*Pi, v=0..Pi);
evalf(%);

4*(1/24 - exp(-64)/24)*Pi

0.5235987756

To obtain a numerical answer you must change the variables to integrate in a cube.
Of course the spherical one is the best here because the new integrand is continuous.
Or, use piecewise:

JJ:=Int(exp(-(sqrt(4*x^2+4*y^2+4*z^2)^3))*piecewise(x^2+y^2+z^2<4,1,0), x=-2..2,y=-2..2,z=-2..2):
evalf(JJ);

0.5235987756

Salut Viorel, bine ai venit!

If you are interested in an approximate solution based on series, here it is:

restart;
ode:= diff(g(r),r$2)- r/R*g(r)=0;
ic:=g(2*R)=0, D[1](g)(0)=R;
dsolve([ode,ic],g(r)) assuming R>0;
G:=rhs(%);
Order:=6;
simplify(series(G, r=0));
convert(%, polynom):
g1:=convert(series(%, R=0), polynom);  
R:=1/3;  #############
dsol:=dsolve([ode,ic],g(r), numeric);
plots:-odeplot(dsol, r=0..2*R);
plot(g1, r=0..2*R);

So, for small R and r, an approximation is

RealDomains is not very reliable but it works in this case (Maple 2019).

eq:=(-2*cos(x)^2+2*sin(x+(1/4)*Pi)^2-1)/sqrt(-x^2+4*x) = 0:
RealDomain:-solve({-x^2+4*x > 0, eq}, {x}, allsolutions);
RealDomain:-solve(eq, {x}, allsolutions);

If you want a compact solution:

9*x-(1/3)*Pi = (-1)^k*(7*x-(1/3)*Pi) + k*Pi:
solve(%,x); # k in Z

I think that in a package where almost everything is a module, such a COMMAND could exist only if it is provided by the package itself.
But we can use:

X := ''RandomVariable''(Normal(mu, sigma)):
Y := eval(X, mu=1);

something like:

plots:-pointplot([x1,y1], color=[seq(COLOR(HUE,(t+1)/2),t=z1)]);

It is not that difficult but unfortunately Maple 2019 gives 0 for it.

It can be computed using residues. This must be done by hand because Maple cannot compute residues for essential singularities.

But it can compute the resulting sum, and the result is   2*Pi*BesselJ(m, b)*exp(I/2*m*Pi)

Download int_iters.mw

On my Windows 7&10, the Maple commands work as expected:

!cmd
                               -1

!calc # ==> opens a console and launch calc.exe (the Windows calculator)
                               0
system(calc);
                               0
system("calc.exe");
                               0
!"cmd /c dir";
                               0
!"cmd /c calc";
ssystem("dir"); # ==>  contents of the current directory (Maple)

restart;
vars := [x,y,z]:
m1,m2,m3 := (a-b)*x*y*z^2, c*x^3, 3*b*x*y*z^2:
EqMon := (m1,m2,vars::list(name)) -> not depends(m1/m2, vars):
EqMon(m1,m2,vars), EqMon(m2,m3,vars), EqMon(m1,m3,vars);
#                       false, false, true

SOL:= u(x, t) = limit(op(1,rhs(sol)),n=0) + subs('n=0..infinity' = 'n=1..infinity', rhs(sol)):
plot3d(subs(infinity=100, rhs(SOL)), x=0..3, t=0..2);

You should inspect visually the integrals.
It's easy to see that the integral wrt phi can be computed symbolically (and actually ==> some integrals are 0).
Then compute numerically wrt theta and r. It will be much faster.

Also, for other 3 dimensional int, try a fast (low precision)  e.g.
Int( ...,  epsilon=1e-2, method=_CubaDivonne)