'maximize'(f) inside the plot command  computes    max {f: t in R}   for a given m  in the interval 1..10.

In Maple this can be expressed:

F:= m -> maximize( cos(2*t/m) + cos(2*(t+5)/m), t );
plot(F, 1..10);

Note that it's not difficult to see that  F(m) = 2*abs(cos(5/m)).
So, the plot can be obtained directly:

plot( 2*abs(cos(5/m)), m=1..10 );

I think that the random behavior is due to the following fact.
For some reason, the order of the terms in a symbolic sum is not deterministic (in some situations); so, because the addition of floats is not associative the result will depend on this order.

add(sin(k),k=1..15);
   sin(1)+sin(2)+sin(3)+sin(4)+sin(5)+sin(6)+sin(7)+sin(8)+sin(9)+sin(10)+sin(11)+sin(12)+sin(13)+sin(14)+sin(15)

add(sin(k),k=1..25);
  sin(19)+sin(20)+sin(21)+sin(22)+sin(23)+sin(24)+sin(25)+sin(2)+sin(3)+sin(4)+sin(5)+sin(6)
  +sin(7)+sin(8)+sin(9)+sin(10)+sin(11)+sin(12)+sin(13)+sin(14)+sin(15)+sin(16)+sin(17)+sin(18)+sin(1)

Note that the order is NOT always random; after a restart it happend that the 1..25 sum became ordered!

J:=proc(n::posint,k::posint)  # numbering 0..n-1
local N;
if n=1 then return 0 fi;
if n<k then return (J(n-1,k)+k) mod n fi;
N:= n-floor(n/k);
floor(k*((J(N,k)-(n mod k)) mod N)/(k-1))
end:

J(5,3);
      3

hterms:=proc(P, n::nonnegint, v::set(name):=indets(P))
  select(t->degree(t,v)=n, expand(P))
end proc:
P := (a^2 + 2)*(b^2 + 2)*(c^2 + 2):
hterms(P,4);

If you think that LambertW is to blame, it is very easy to get rid of it replacing in your relation

LambertW(x)=y   with  x = y*exp(y).

But probably the problem is that in your expression the rounding errors are too big. Note that if Digits is increased to 30 the plot is almost empty, so this could be the real situation.

Dr:=proc(P::depends(list(And(ratpoly(integer,Variables),Not(integer)))),Variables::list(symbol),DEvar::symbol,DEsuffix::string)

PS. Note that ratpoly means rational function (such as (x+1)/(x*y+2)), not polynomial.

It is enough to set

_EnvFormal:=true;

Note that without this, even sum(p^k, k=0..infinity) is not computed. Of course one may use assume on parameters (in order to have convergence) but there are too many of them and I understand that your series are formal.

convert(a, cot);

I'd recommend to start with these much simpler integrals and see what happens:

int(sqrt(x^4+1),x);
int(sqrt(x^4+1+x),x);
int(sqrt(x^4+1+sqrt(x^2+1)),x);

restart;
with(LinearAlgebra):
t:=time():
n:=1000:
A:=RandomMatrix(n,n,generator=rand(-10.0..10.0),shape=symmetric,datatype=float[8]):
V,Q:=Eigenvectors(A):
MV:=DiagonalMatrix(sqrt~(Vector(V,datatype=complex))):
B:=Q.MV.Q^*:
`||A-B^2||`=Norm(A-B^2);
'time'=time()-t;

                                           
             ||A-B^2|| = 4.49369731341387095 * 10^(-11)  
                         time = 26.021

The integral `f*g` is expressed with the Dirac "function" and some problems appear.
Just convert f to Heaviside and it will work.
You may insert after your f:=...

f:=unapply( convert(f(x),Heaviside), x );

(Note that `f*g`(x) will be much simpler!)

Now OK :

TestCylinderInscritSphere-ok.mw

B:=copy(A);
Or,
B:=Matrix(A);

It is easy to see that this matrix was specially cooked to have no sqrt.

Its Jordan form is:

JordanForm(A);

but the first 2x2 block  [[0,1], [0,0]] does not have a square root, so A neither.

Edit.  Unfortunately I made a stupid (but hard to detect) mistake. The fact that a Jordan block does not have a sqrt does not imply that that the  whole Jordan form does not have sqrt; only the converse is true.

In this case, sqrt(A) exists (not unique of course) and can be computed. It is:

B:=Matrix(9, 9, [[2203/369+(41/135)*sqrt(5), 8665/1476+(14/45)*sqrt(5), 1145/246+(43/135)*sqrt(5), -1687/492+(44/135)*sqrt(5), 28/41+(1/3)*sqrt(5), -3979/1476+(46/135)*sqrt(5), -1315/738+(47/135)*sqrt(5), 125/492+(16/45)*sqrt(5), -3515/369+(49/135)*sqrt(5)], [-11015/1476+(41/135)*sqrt(5), -102365/5904+(14/45)*sqrt(5), 917/984+(43/135)*sqrt(5), 7451/1968+(44/135)*sqrt(5), -99/164+(1/3)*sqrt(5), 61223/5904+(46/135)*sqrt(5), -14089/2952+(47/135)*sqrt(5), 6263/1968+(16/45)*sqrt(5), 17575/1476+(49/135)*sqrt(5)], [4753/1476+(41/135)*sqrt(5), 102187/5904+(14/45)*sqrt(5), -5875/984+(43/135)*sqrt(5), -5149/1968+(44/135)*sqrt(5), -11/164+(1/3)*sqrt(5), -70897/5904+(46/135)*sqrt(5), 28319/2952+(47/135)*sqrt(5), -5281/1968+(16/45)*sqrt(5), -10001/1476+(49/135)*sqrt(5)], [-9539/1476+(41/135)*sqrt(5), -97937/5904+(14/45)*sqrt(5), 1409/984+(43/135)*sqrt(5), 7943/1968+(44/135)*sqrt(5), -99/164+(1/3)*sqrt(5), 59747/5904+(46/135)*sqrt(5), -15565/2952+(47/135)*sqrt(5), 4787/1968+(16/45)*sqrt(5), 16099/1476+(49/135)*sqrt(5)], [7705/1476+(41/135)*sqrt(5), 111043/5904+(14/45)*sqrt(5), -4891/984+(43/135)*sqrt(5), -4165/1968+(44/135)*sqrt(5), -11/164+(1/3)*sqrt(5), -73849/5904+(46/135)*sqrt(5), 25367/2952+(47/135)*sqrt(5), -8233/1968+(16/45)*sqrt(5), -12953/1476+(49/135)*sqrt(5)], [-5111/1476+(41/135)*sqrt(5), -72845/5904+(14/45)*sqrt(5), 1901/984+(43/135)*sqrt(5), 5483/1968+(44/135)*sqrt(5), 65/164+(1/3)*sqrt(5), 55319/5904+(46/135)*sqrt(5), -19993/2952+(47/135)*sqrt(5), 359/1968+(16/45)*sqrt(5), 11671/1476+(49/135)*sqrt(5)], [10657/1476+(41/135)*sqrt(5), 113995/5904+(14/45)*sqrt(5), -2923/984+(43/135)*sqrt(5), -7117/1968+(44/135)*sqrt(5), 153/164+(1/3)*sqrt(5), -64993/5904+(46/135)*sqrt(5), 16511/2952+(47/135)*sqrt(5), -9217/1968+(16/45)*sqrt(5), -15905/1476+(49/135)*sqrt(5)], [-11015/1476+(41/135)*sqrt(5), -137789/5904+(14/45)*sqrt(5), 4853/984+(43/135)*sqrt(5), 7451/1968+(44/135)*sqrt(5), -99/164+(1/3)*sqrt(5), 84839/5904+(46/135)*sqrt(5), -25897/2952+(47/135)*sqrt(5), 10199/1968+(16/45)*sqrt(5), 17575/1476+(49/135)*sqrt(5)], [4753/1476+(41/135)*sqrt(5), 49051/5904+(14/45)*sqrt(5), 29/984+(43/135)*sqrt(5), -5149/1968+(44/135)*sqrt(5), -11/164+(1/3)*sqrt(5), -35473/5904+(46/135)*sqrt(5), 10607/2952+(47/135)*sqrt(5), 623/1968+(16/45)*sqrt(5), -10001/1476+(49/135)*sqrt(5)]]);

# Proof:
Norm(A - B^2);
        0

delta is supposed to be in 0 .. 1

Then max(min(...)) = 1 - delta.

==> lim = 1 - delta

[You must use maths, not Maple for this].