Straight from the Maple worksheet window:

Help -> On the Web -> Application Center

and download the applications that interest you. 

Unfortunately for the problem to collect  rhs(aa)  for  omega^(-sigma+1)*L   neither  frontend  nor  freeze  method does not work. Here 2 more complicated ways to solve the problem.

First way:

lhs(aa)=applyop(collect, [1,1], frontend(collect, [rhs(aa), omega^(-sigma+1)]), L);

Second way uses  subs  and   collect  commands only:

lhs(aa)=subs(t=omega^(-sigma+1)*L,collect(subs(omega^(-sigma+1)=t/L, rhs(aa)), t));

The solution to the second problem:

lhs(aa)=applyop(factor, [1,1], frontend(collect, [rhs(aa), omega^(-sigma+1)]));

or simply

simplify(aa);

for full factorization.

restart;

r:=6:
plot3d([r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)],  phi=0..2*Pi,  theta=0..Pi);

Addition: The shortest way of plotting a sphere is

restart;
plots[display](plottools[sphere]([0,0,0], 6));

Small circles at the ends of the segments should be plotted separately:

f:=x->piecewise(type((x+2)/4,integer)=false,x-floor((x+2)/4)*4, undefined):  # Setting for your function
f(2.5), f(-2), f(2), f(6);
A:=plot(f, -10..10.1, scaling=constrained, discont):
B:=plot([seq([-10+4*k,-2],k=0..4), seq([-6+4*k,2],k=0..4) ],style=point, symbol=circle):
plots[display](A, B);

f(x,y,z)=1  is a surface in the 3D-space. It can be plotted by  plots[implicitplot3d]  command.

Example:

plots[implicitplot3d](x^2+y^2-z^2=1, x=-3..3, y=-3..3, z=-3..3,  numpoints=8000);

Use  plottools[getdata]  and  plots[display]  with  insequence  option.

Example:

A:=plots[implicitplot](x^5+y^5-x*y=1, x=-3..3,y=-3..3, gridrefine=2):
B:=plottools[getdata](A):
C:=convert(B[3],listlist):
plots[display](seq(plot(C[1..n]), n=1..nops(C), 3), insequence);

You can use  seq  command:

DataName:=Array(["Data1", "Data2", "Data3"]):
DataSet:=[seq(i=DataName[i], i=1..numelems(DataName))];

                      DataSet := [1 = "Data1", 2 = "Data2", 3 = "Data3"]

Edited.

In this line 

fplt[j] := plots[odeplot](sol2, sol1, [[eta, diff(f(eta), eta)], [eta, diff(f(eta), eta)]], color = K1[j], axes = boxed);

should be one thing  only  sol2  or  sol1

Combine these graphs on one plot you can by  plots[display]  command.

By this way you can build angles on the plane anywhere and any values:

restart;
with(plottools): with(plots):
A:=line([0.5,0.5],[1.5,0.5], color=red, thickness= 3):
B:=rotate(A, 9*Pi/4, [0.5,0.5]):
C:=plot([0.5+(0.1+0.015*t)*cos(t),0.5+(0.1+0.015*t)*sin(t), t=0..9*Pi/4], color=black):
P:=curve([[0.655,0.6],[0.647,0.647],[0.68,0.625]]):
T:=textplot([0.45,0.45, 9*Pi/4], color=blue):
display(A, B, C, P, T, scaling=constrained, view=[0..1.7,0..1.5]);

Edited. Arrow added at the end of the arc, which shows the direction of the rotation.

If you want to immediately build a few graphs for specific values of the parameter m, and of specific colors, you can do so:

M:=[0.5, 2, 5, 8, 10]:  Color:=[red,green,blue,gold,black]:
plot([seq(-arctan(2*m*x/(1-x^2)), m=M)], x=0..10, color=Color,  numpoints=1000, discont, scaling=constrained);

All graphs have finite discontinuities at  x = 1  ( limit(f(x), x=1, left)=-Pi/2  and  limit(f(x), x=1,right)=Pi/2 )

Carl, you have incorrectly entered the expression for the function.

Pretty good approximation gives the model function with  arctan :

x_val:=<250,300,350,397,451,497,547,593,647,691,745,788,840,897>:

y_val:=<0,0.5,2,6.3,23.2,48.7,71.2,83.4,90.1,92.8,94.7,95.7,96.9,97.8>:

A:=plot(x_val, y_val, style=point, symbol=solidcircle, color=blue, symbolsize=12):  # Plot of initial data

f:=Statistics[NonlinearFit](a+b*arctan(c*x+d),x_val, y_val, x, initialvalues=[a = 50, b = 33, c = 0.02, d=-12]);  # Model function

B:=plot(f, x=x_val[1]..x_val[-1], color=red): # Plot of the model function

plots[display](A, B);

Addition: If the values of the y-variable are the percentages, then we can simply scale these values. The graph will have the same form:

x_val:=<250,300,350,397,451,497,547,593,647,691,745,788,840,897>:
y_val:=1/100*<0,0.5,2,6.3,23.2,48.7,71.2,83.4,90.1,92.8,94.7,95.7,96.9,97.8>:
A:=plot(x_val, y_val, style=point, symbol=solidcircle, color=blue, symbolsize=12):  # Plot of initial data
f:=Statistics[NonlinearFit](a+b*arctan(c*x+d),x_val, y_val, x, initialvalues=[a = 50/100, b = 33/100, c = 0.02, d=-12]);  # Model function​​​​​​​
B:=plot(f, x=x_val[1]..x_val[-1], color=red): # Plot of the model function
plots[display](A, B);

A more detailed investigation shows that  f  has not the values  -1 or 1 anywhere.  Therefore, the range  of this function is the open interval  -1 .. 1 :

f := sqrt(x^2+x+1)-sqrt(x^2-x+1):

solve(diff(f,x)>0);  # The function is strictly increasing on the whole axis  x

limit(f, x=-infinity);

limit(f, x=infinity);

minimize(f, location); 

maximize(f, location);

plot(f, x=-10..10);

Joe, your first method is very witty, vote up! Here is another simple method using ListTools  package tools. It does exactly what the OP wants to:

L := op~([g[ ]]);
S:={L[ ]};
for x in S do print(cat(x,`  appear  `, nops([ListTools[SearchAll](x, L)]))); end do;
 

A:=2*Pi^2*(-1+delta)/r^(1+delta);
-evalindets(A, 'And(`+`, satisfies(f->sign(op(1,f))=-1))', t->-t);

Tom, your procedure will not work if one of the ends of the range is equal to infinity. Here's my version:

ConvertTo_range:=proc(S::{set,list})

local T, T1;

T:=map(t->[op(t)],S);

if nops(T)=2 then

T1:=seq(op(remove(s->type(s,name),T[i])), i=1..2);

return min(T1)..max(T1) else

if type(op(T)[1],name) then return -infinity..op(T)[2] else op(T)[1]..infinity fi; fi;

end proc:

Example of use:

L:=[{ arccos(3/4) < theta},

{arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2)))+Pi <= theta, theta < 2*Pi-arccos(3/4)},

{theta <= 3*Pi-arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2))), 2*Pi+arccos(3/4) < theta}]:

ConvertTo_range~(L);

evalf(%);