restart;
r := 1.9;
i := floor(r);

You can use  Student:-Calculus1:-Roots  command to find all the real roots in the specified interval. It does not need  maxsols  option and can work symbolically.

Examples:

Student:-Calculus1:-Roots(eq, beta=-10..10);
Student:-Calculus1:-Roots(sin(x), x=-10..10);
 

Of course in the real domain it is better to write  int(tan(x), x) = -ln(abs(cos(x)) . The same for the integrals  int(1/x, x) , int(1/(x^2-1), x)  and so on. But the result of Maple is also correct, since it works by default in the complex domain and the results may differ by a complex constant. See

simplify([ln(-2), ln(-5)]);
                                           [ln(2)+I*Pi, ln(5)+I*Pi]

Since, as you can see, the results have the same complex constant, the final result is correct when calculating certain integrals using the Newton – Leibniz formula:

A:=int(tan(x), x);
B:=int(tan(x), x=2*Pi/3..3*Pi/4);
C:=``(expand(eval(A, x=3*Pi/4))) - ``(eval(A, x=2*Pi/3));
expand(C);
              
                                               

The first result  B  ( cos(x)<0  in this range  2*Pi/3..3*Pi/4 ) is calculated by Maple, and the second one  C  is calculated directly by the Newton-Leibniz formula. The results are the same and correct.

Probably a more complete answer why Maple is designed in this way, the developers of  int  procedure can give you.

Unfortunately Maple hangs when trying to compute it symbolically. However, for practical applications, a numerical solution is sufficient. The procedure  V  numerically finds the volume of intersection of two spherical caps (your picture below).

restart;
V:=proc(r,h,h__a)
local r0;
Digits:=20;
r0:=sqrt(2*r*h-h^2);  # This is the base radius of the top cap
if not (is(h>0) and is(h<=r) and is(h__a>r-r0) and is(h__a<=r)
) then error "Invalid arguments" fi;
evalf(Int(1, [z=r-h..sqrt(r^2-x^2-y^2),y=-sqrt(r0^2-x^2)..sqrt(r0^2-x^2),x=r-h__a..r0]));
end:

Examples of use:

V(1, 0.8, 0.6);
V(1, 0.1, 0.5);
V(1, 1, 1);  # This is a quarter ball radius 1
evalf(Pi/3);  # Check

Edit.

I simplified your code by removing unnecessary indexes and replacing the for-loop with  seq  command. You can check that there are no real solutions for thetabn > 43*Pi/100 :

restart;
v := 145000;
thetavn := (1/6)*Pi;
omegac := 0.1;
s := v*cos(thetavn)*(cos(2*thetabn)*tan(thetabn)+sin(2*thetabn)*sin(omegac*t)/omegac);
 Data:=[seq([thetabn, solve(s = 0, t)], thetabn=[seq(Pi*k/100, k=1..43)])];
plot(Data);

You can define expr using the  piecewise  command:

expr:=piecewise(b<>a^3, (a-b)/(b-a^3), undefined); 

Examples of use:
eval(expr, [a=2,b=8]);
eval(expr, [a=2,b=7]);
 

Unfortunately Maple can simplify this only for a specific  L  and a specific  f :

restart; with(IntegrationTools):
L:=3:
int(f(x), x = 0 .. L*Ts)-Split(int(f(x), x = 0 .. L*Ts), [seq(i*Ts, i = 1 .. L-1)]);
eval(%, f=sin);
       

Here it works:

is(5^(x-1)=5^x/5);
                                           true

First we solve this system as a system of equations, leaving only one inequality  p[1, 2] > 0 . We see that the solution depends on only one positive parameter p[1, 2] (the straight line from the solutions in 3D). So easy to get a general solution:

S:={0 < p[1, 2], p[1, 1] < 2*p[2, 2]+(3/2)*p[1, 2], p[1, 2]^2/p[2, 2] < p[1, 1], (2/3)*p[1, 2] < p[2, 2]};
solve({0 < p[1, 2], p[1, 1] <= 2*p[2, 2]+(3/2)*p[1, 2], p[1, 2]^2/p[2, 2] = p[1, 1], (2/3)*p[1, 2] = p[2, 2]});
{p[1, 1] = 3*p[1, 2]*(1/2), p[2, 2] > 2*p[1, 2]*(1/3), 0 < p[1, 2]};  # This is the general solution
eval(S, {p[1,2]=6, p[1,1]=9, p[2,2]=5});  # Check
              

Download code_new.mw

Or replace the last line with:

work := simplify(work, {G*M/R^2 = g});
                                                                   work := m*h*g

It works properly in Maple 2018.2 . I believe this is the issue of Maple's version.

Addition. Just try

plot3d(sin(x)*cos(y), x = 0 .. 4*Pi, y = 0 .. 4*Pi, view = [default, default, -3 .. 3],  color=x);
 

If you want to find the maximum coefficient of a polynomial (PSI  is a polynom of X) and its position i (the coefficient in front of  X^i ), then it’s better not to use a looping search, but  max[index]  command. Here is the solution to your example:

restart;
G := unapply((x-y+1)*X/binomial(x, y)+1-(x-y+1)/binomial(x, y), x, y); 
PSI := G(3, 2)^10*G(4, 2)^4*G(4, 3)^2*G(5, 2)^2*G(5, 4)*G(5, 3)*G(6, 2)^4*G(7, 3)*G(8, 2)*G(9, 2)*G(9, 3)*G(9, 4)*G(10, 2)^2*G(11, 2)*G(11, 3)*G(12, 3)*G(13, 3)*G(15, 2)^2;
PSI:=evalf(expand(PSI));
L:=[seq(coeff(PSI,X,i), i=0..degree(PSI))];
max[index](L);
%-1, L[%];  # The final result
                                                  14,  0.1499364076

I removed some of the lines that do not affect anything, and also made small changes in the last lines to get graphics of acceptable quality:

restart:
with(plots):

Omega:=5*Pi:
gamma1:=8*Pi:
gamma2:=0.002*Pi:
x1:=100*Pi:
omega2:=200*Pi:
gamma0:=0.2*Pi:
G:=20*Pi:

lambda1:=(1/(2*Pi))^2*(G*Omega*gamma0/(2*(0.25*gamma0^2+Delta^2))):
lambda2:=(1/(2*Pi))^2*(-G*Omega*Delta/((0.25*gamma0^2+Delta^2))):
lambda3:=(1/(2*Pi))^2*gamma1+lambda1:
lambda4:=(1/(2*Pi))^2*(0.2*omega2-G^2*Delta/(0.25*gamma0^2+Delta^2)):
lambda5:=(1/(2*Pi))^2*(2*x1^2*omega2/((omega2^2+gamma2^2))):
f:=epsilon-(lambda1+sqrt(-lambda2^2+Y*(lambda3^2+(lambda4-lambda5*Y)^2)))^2:

implicitplot3d(f, epsilon=0..20,Y = 0 .. 0.4,Delta=-25*Pi..25*Pi, style=surface, color=khaki, labels=[E1,Y,Delta],tickmarks=[3,3,3], numpoints=500000);

contourplot(solve(f,epsilon),Delta =-25*Pi .. 25*Pi,Y=0..1,contours=[4,8,12,16,20],axes=boxed,thickness=2,color=black,font=[1,1,18],tickmarks=[5, 5],linestyle=1, numpoints=50000);

The procedure called  GM  solves your problem:

GM:=proc(L,B)
local coefff, T;
T:=indets(B);
coefff:=proc(P,T,t) 
local C,H,i,k:
C:=[coeffs(P,T,'h')]: H:=[h]: k:=0:
for i from 1 to nops(H) do
if H[i]=t then k:=C[i] fi:
od:
k;
end proc;
map(p->Matrix(nops(p),nops(B), (i,j)->coefff(p[i],T,B[j])), L);
end proc:

#  Example of use:

L:=[[a*x^2+b*x*y-1, -(a*b-b)*x*y/a-z^2+(a-1)/a, -a*c*z^2/(b*(a-1))+(b+c)/b], [a*x^2+b*x*y-1, -(a*b-b)*x*y/a-z^2+(a-1)/a, 1],[a*x^2+b*x*y-1, -z^2+(a-1)/a, c*x*y+1], [b*x*y-1, -x^2-z^2, (b+c)/b], [-1, -x^2-z^2]]:
B:=[x^2,x*y,z^2,1]:
GM(L, B);