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)

s4:=series(ex, epsilon, 4): # ex is your expression
coeff(s4, epsilon, 2); # coeff of epsilon^2

A := Matrix(2,3,()->randpoly(t,degree=1));
int~(A,t=-1..1);

That is what subs must do:

subs(a=b, {c=d}, {e=f}, {g=h} ) ==> successive substitutions  a=b, c=d, e=f  in  {g=h}.
In your case sol := {c=d}, {e=f}, {g=h}

is_solution_trivial:= proc(pde, sol::{`=`, set(`=`),`function`})
try
  evalb(eval(`if`(sol::set, sol, {sol}), depvars(pde)=~ 0) = {0=0})
  catch : false   # or catch "numeric exception: division by zero"
end try
end proc;