Use initial conditions and  numeric  option. You can get the values of a function at individual points and plot its graph.

Example:

Sol:=dsolve({-1.699678499+diff(z(t), t, t)-3.813322790*z(t)+.6966019659*z(t)^2+0.1548004368e-1/z(t)^2 = 0, z(0)=1, D(z)(0)=0}, z(t), numeric);
Sol(1);
plots:-odeplot(Sol, [t,z(t)], t=0..10, view=0..10);

Obviously this is bug. Even if you use  RealDomain:-solve  command, the error persists :

RealDomain:-solve(-sqrt(k*x) = sqrt(k*x+k-1), {k,x});

                               {k = 1, x = x}

z:=x+I*y:
evalc(abs(z+3-2*I) + abs(z-3-8*I)) = 6*sqrt(2);
sqrt((3+x)^2+(y-2)^2)+sqrt((-3+x)^2+(y-8)^2) = 6*sqrt(2)

Geometrically, this equation defines on the plane  (x, y)  the set of such points, the sum of the distances from each of which to fixed points  A(-3, 2)  and  B(3, 8)  is equal to  6*sqrt(2) . But the distance between these points is the same number  6*sqrt(2). Therefore this set is the segment joining points  A  and  B : .  Therefore, the minimum and maximum is reached inside or at the ends of this segment:

A:=[-3,2]: B:=[3,8]:
C:=(1-t)*~A+t*~B:  # Arbitrary point of the segment  AB
simplify([(minimize,maximize)(sqrt(C[1]^2+C[2]^2), t=0..1)]);

                        [5*sqrt(2)*(1/2), sqrt(73)]

Edit.


                                     

Here is another very simple way, which is even slightly more effective than  select :

CodeTools:-Usage(seq(`if`(isprime(n),n,NULL),n=10^5..10^6)):
memory used=231.00MiB, alloc change=-3.43MiB, cpu time=4.67s, real time=4.72s, gc time=296.88ms

CodeTools:-Usage(select(isprime, [$10^5..10^6])):
memory used=268.77MiB, alloc change=3.43MiB, cpu time=5.31s, real time=5.31s, gc time=687.50ms


 

I don't know such the function, but the simple procedure  PrimesInRange  solves the problem:

PrimesInRange:=proc(a, b)
local p, n;
p[1]:=nextprime(a-1); 
for n from 1 while p[n]<=b do
p[n+1]:= nextprime(p[n]) 
od:
seq(p[k], k=1..n-1);
end proc:

Example of use:

PrimesInRange(1, 100);
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97

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