I assumed that  S>-Pi/2  and  S<Pi/2 :

f:= [((2*a*sin(S)*cos(S)^(2)))/(1-sin(S)^(3))<1,-Pi/2<S,S>Pi/2,a>1];
A:=plots:-implicitplot([2*a*sin(S)*cos(S)^2/(1-sin(S)^3) = 1, a=1], S=-Pi/2..Pi/2, a=-5..10, color=[red, blue], thickness=2, gridrefine=3):
B:=plots:-inequal({f[1], f[4]}, S=-Pi/2..Pi/2, a=-5..10, optionsfeasible = [color = "LightGreen"]):
plots:-display(A,B);

The set of solutions corresponds to the green region, which is unbounded from above.

Of course, this system is easy to solve purely analytically, if we consider the cases  S<0, S=0, S>0 .

The results obtained with the help of  extrema  command require additional justification, since among them there can be also saddle points in which neither a maximum nor a minimum is reached, as in the example:

extrema(x^2-y^2, { }, {x,y}, 's'); 
s;  # This is a saddle point
                                {0}
                        {{x = 0, y = 0}}

Also, extrema  can not return results at all, if the differentiability is violated at extreme points, as in the example below ((0,0) is the point of the global minimum):
extrema(sqrt(x^2+y^2), { });  # NULL is returned

As an alternative to extrema you can use Optimization:-Minimize  and  Optimization:-Maximize commands:

z:=a+b*I:                                                     
Digits:=30:
fnormal(Optimization:-Minimize(evalc(abs(z)), evalc({abs(z+1)+4*abs(z-1)=25})), 20);
fnormal(Optimization:-Maximize(evalc(abs(z)), evalc({abs(z+1)+4*abs(z-1)=25})), 20);

 [4.4000000000000000000, [a = -4.4000000000000000000, b = -0.]]
  [5.6000000000000000000, [a = 5.6000000000000000000, b = 0.]]

Of course, these commands are not ideal and in general they can return not global, but local extremes.
 

Use  PolyhedralSets  package (was introduced in Maple 2015). 

The above example:

with(PolyhedralSets):
V:=[[2,1,0], [-2,1,0], [-2,-1,0], [2,-1,0], [2,1,0], [1,0,2], [2,-1,0], [2,1,0], [-2,1,0], [-1,0,2], [1,0,2],  [-2,1,0], [-1,0,2], [-2,-1,0], [2,-1,0], [-2,-1,0], [-1,0,2], [1,0,2]]:  # The list of all the vertices
ps := PolyhedralSet(V);
Plot(ps);
Faces(ps);
Vertices(ps);
Edges(ps);

If you have an older Maple, then write. In this case, a procedure is necessary to solve your problems.

It seems  combine  command makes it:

combine(16*sin(x)^2*cos(y)^3);

         -cos(2*x-3*y)-cos(2*x+3*y)-3*cos(2*x-y)-3*cos(2*x+y)+2*cos(3*y)+6*cos(y)

If all faces of a polyhedron are given as lists of lists the coordinates of its vertices, then it is easy to construct this polyhedron using the commands  plottools:-polygon  and  plots:-display . 

Example:

Faces:=[[2,1,0],[-2,1,0],[-2,-1,0],[2,-1,0]], [[2,1,0],[1,0,2],[2,-1,0]], [[2,1,0],[-2,1,0],[-1,0,2],[1,0,2]], [[-2,1,0],[-1,0,2],[-2,-1,0]], [[2,-1,0],[-2,-1,0],[-1,0,2],[1,0,2]]:
S:=seq(plottools:-polygon(Faces[i], color=khaki), i=1..5):
plots:-display(S, scaling=constrained);

Unfortunately I do not know how to make it easier:

restart;
x__1:=y^(1/(a+b))/(w__1*b/(a*w__2))^(b/(a+b));
expand(x__1) assuming positive;
B:=combine(%);
f:=c->applyrule(x::anything^(-y::anything)*z::anything^y::anything=(z/x)^y, c);
f(B);
applyop(expand, 1, %) assuming a>0, b>0, w__1>0, w__2>0;

For  x2  everything is the similar.

The variable  v  determines the width of the Mobius strip. If  v = -1  then we get half the edge, and  v = 1  gives the other half. For details see wiki  https://en.wikipedia.org/wiki/M%C3%B6bius_strip
 

L:=[(1+v/2*cos(u/2))*cos(u), (1+v/2*cos(u/2))*sin(u), v/2*sin(u/2)]:
S:=plot3d(L,  u = 0 .. 2*Pi, v=-1..1): 
C1:=plots:-spacecurve(eval(L, v=-1), u = 0 .. 2*Pi, color=red, thickness=3):
C2:=plots:-spacecurve(eval(L, v=1), u = 0 .. 2*Pi, color=red, thickness=3):
plots:-display(C1, C2, S, title = "Möbius strip", titlefont = [Courier, bold, 14], axes=normal, scaling = constrained, orientation=[-25,65]); 

Addition:  animation of the bypass of the edge of the Mobius strip:

restart;
L:=[(1+v/2*cos(u/2))*cos(u), (1+v/2*cos(u/2))*sin(u), v/2*sin(u/2)]:
S:=plot3d(L, u = 0 .. 2*Pi, v=-1..1, style=surface): 
plots:-animate(plots:-spacecurve,[eval(L,v=1.02),u = 0 .. a, color=red, thickness=4, axes=normal, scaling = constrained, orientation=[-25,65]], a=0..4*Pi, frames=60, background=S);

collect command does not work because only one member in esp1 contains  a^2. You can carry out your decomposition by hand:

esp1 := -a^4+a^2*c^2:
a^2*simplify(esp1/a^2);

                                 a^2*(-a^2+c^2)

The main difference is that  assuming  imposes a restriction only when this command is executed:

Examples:

restart;
simplify(sqrt((x-2)^2)) assuming x<=2;
about(x);

                                    2 - x
            x:  nothing known about this object

assume(x<=2):
simplify(sqrt((x-2)^2));
about(x);
                                    2 - x
Originally x, renamed x~:
  is assumed to be: RealRange(-infinity,2)

 

For more information, see Help on these commands.

Add this line at the end of your code:

evalm(Ltot);

Since your expression is a square trinomial with respect to x, it is sufficient to solve the system  {f(a,b)>0, discrim(Delta,x)<0} 

Example:

restart;
Delta:=f(a,b)*x^2+g(a,b)*x+h(a,b):
Sys:={f(a,b)>0, discrim(Delta,x)<0};
f(a,b):=1-a:  g(a,b):=a+b:  h(a,b):=a-b:
solve(Sys);

See help on  ?procname(recursion)

We need to help Maple a bit:

f:=(5*x^2+3*(x+exp(x))^(1/3)+exp(x)*(2*x^2+3*x))/x/(x+exp(x))^(1/3);
IntegrationTools:-Expand(int(f, x));
simplify(%);  # The result
simplify(diff(%, x));  # Checking

Addition. In fact, the second integral  int(1/(sqrt(x+exp(x))), x)  is more hard. Mathematica fails with it also.

For  fsolve  command, the number of equations must equal the number of unknowns. I took  V=1 :

restart;
omega1 := 1.562:
omega2 := 2.449:
omega3 := 3.325:
y1 := c1*sin(omega1*t+phi1)+c2*sin(omega2*t+phi2)+c3*sin(omega3*t+phi3): 
y2 := .1019*c1*sin(omega1*t+phi1)+.75*c2*sin(omega2*t+phi2)+.4608*c3*sin(omega3*t+phi3):
y3 := .407*c1*sin(omega1*t+phi1)+(0*c2)*sin(omega2*t+phi2)+1.844*c3*sin(omega3*t+phi3): 
eq1 := subs(t = 0, y1) = 0:
eq2 := subs(t = 0, y2) = 0: 
eq3 := subs(t = 0, y3) = 0: 
eq4 := subs(t = 0, diff(y1, t)) = V:
eq5 := subs(t = 0, diff(y2, t)) = 0:
eq6 := subs(t = 0, diff(y3, t)) = 0:
V:=1: 
eqs := {eq1, eq2, eq3, eq4, eq5, eq6}: 
vars := {c1, c2, c3, phi1, phi2, phi3}: 
fsolve(eqs, vars);

    {c1 = -0.8218023979, c2 = 0.0001356814559, c3 = 0.08520994334, phi1 = -9.424777961, phi2 =    3.141592654, phi3 = -3.141592654}

You can save V variable without assigning any value to it. In this case, you must use  solve  command and you will get several solutions:

restart;
omega1 := 1.562:
omega2 := 2.449:
omega3 := 3.325:
y1 := c1*sin(omega1*t+phi1)+c2*sin(omega2*t+phi2)+c3*sin(omega3*t+phi3): 
y2 := .1019*c1*sin(omega1*t+phi1)+.75*c2*sin(omega2*t+phi2)+.4608*c3*sin(omega3*t+phi3):
y3 := .407*c1*sin(omega1*t+phi1)+(0*c2)*sin(omega2*t+phi2)+1.844*c3*sin(omega3*t+phi3): 
eq1 := subs(t = 0, y1) = 0:
eq2 := subs(t = 0, y2) = 0: 
eq3 := subs(t = 0, y3) = 0: 
eq4 := subs(t = 0, diff(y1, t)) = V:
eq5 := subs(t = 0, diff(y2, t)) = 0:
eq6 := subs(t = 0, diff(y3, t)) = 0: 
eqs := {eq1, eq2, eq3, eq4, eq5, eq6}: 
vars := {c1, c2, c3, phi1, phi2, phi3}: 
solve(eqs, vars);

restart;
evalc((2-I*X)^4);
Re(%) assuming X::real;

or

restart;
expand(Re((2-I*X)^4)) assuming X::real;

or

restart;
evalc(Re((2-I*X)^4));

Probably the last method is the shortest  (evalc  assumes that all parameters are real numbers).