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);

Unfortunately Maple 2017 fails with this integral (for Maple 2018 no problems). You can easily calculate the integral in Maple 2017  in polar coordinates by writing it as a double integral:

int(eval(2 - 2*sqrt(x^2+y^2), {x=r*cos(t), y=r*sin(t)})*r, [r=0..1, t=0..2*Pi]); 
                                             2*Pi*(1/3)

S := [x^2 = A1, x*y = A2, z^2 = A3] :
SS := map(t->rhs(t)=lhs(t), S);

You forgot to call  LinearAlgebra  package:

A := `<|>`(`<,>`(-1, -3, -6), `<,>`(3, 5, 6), `<,>`(-3, -3, -4));
v, e := LinearAlgebra:-Eigenvectors(A);