Perhaps for the inexperienced beginners, the following code will be more understandable (on the example of the initial list):

restart;
R:=[4,-1,9,11,-4,2,4,1,-1,2.2,112,44,-134,0.124,34,41,7,97,13.22123]:
n:=nops(R):
for i from 1 to n do
v:=R[i];
if v>-2 and v<3 then vals[i]:=R[i]; Numbs[i]:=i  fi;
od:
vals:=convert(vals,list);
Numbs:=convert(Numbs,list);

                  vals := [-1, 2, 1, -1, 2.2, 0.124]
                    Numbs := [2, 6, 8, 9, 10, 14]
 

The system has a unique solution:

Sys:={seq(x[j-1]-2*x[j]+x[j+1]=0, j=3...98),
-2*x[1]+x[2]=5,
x[1]-2*x[2]+x[3]=-4,
x[98]-2*x[99]+x[100]=-8,
x[99]-2*x[100]=13}:
fsolve(Sys);

{x[1] = -1., x[2] = 3., x[3] = 3., x[4] = 3., x[5] = 3., x[6] = 3., x[7] = 3., x[8] = 3., x[9] = 3., x[10] = 3., x[11] = 3., x[12] = 3., x[13] = 3., x[14] = 3., x[15] = 3., x[16] = 3., x[17] = 3., x[18] = 3., x[19] = 3., x[20] = 3., x[21] = 3., x[22] = 3., x[23] = 3., x[24] = 3., x[25] = 3., x[26] = 3., x[27] = 3., x[28] = 3., x[29] = 3., x[30] = 3., x[31] = 3., x[32] = 3., x[33] = 3., x[34] = 3., x[35] = 3., x[36] = 3., x[37] = 3., x[38] = 3., x[39] = 3., x[40] = 3., x[41] = 3., x[42] = 3., x[43] = 3., x[44] = 3., x[45] = 3., x[46] = 3., x[47] = 3., x[48] = 3., x[49] = 3., x[50] = 3., x[51] = 3., x[52] = 3., x[53] = 3., x[54] = 3., x[55] = 3., x[56] = 3., x[57] = 3., x[58] = 3., x[59] = 3., x[60] = 3., x[61] = 3., x[62] = 3., x[63] = 3., x[64] = 3., x[65] = 3., x[66] = 3., x[67] = 3., x[68] = 3., x[69] = 3., x[70] = 3., x[71] = 3., x[72] = 3., x[73] = 3., x[74] = 3., x[75] = 3., x[76] = 3., x[77] = 3., x[78] = 3., x[79] = 3., x[80] = 3., x[81] = 3., x[82] = 3., x[83] = 3., x[84] = 3., x[85] = 3., x[86] = 3., x[87] = 3., x[88] = 3., x[89] = 3., x[90] = 3., x[91] = 3., x[92] = 3., x[93] = 3., x[94] = 3., x[95] = 3., x[96] = 3., x[97] = 3., x[98] = 3., x[99] = 3., x[100] = -5.}

with(DynamicSystems):
T:=Vector(4, t->t):
A:=Vector(4, t->1+sin(t)):
DiscretePlot(evalf(T), evalf(A), style=stem);

The intersection of two cylindars:

B1:=plots:-spacecurve(1.005*[cos(t),sin(t),cos(t)], t=0..2*Pi, color=blue,thickness=3):
B2:=plots:-spacecurve(1.005*[cos(t),sin(t),-cos(t)], t=0..2*Pi, color=blue,thickness=3):
S1:=plot3d([[cos(t),sin(t),z],[-cos(t),sin(t),z]], t=-Pi/2..Pi/2, z=-cos(t)..cos(t), style=surface, color=green):
S2:=plot3d([[x,cos(t),sin(t)],[x,cos(t),-sin(t)]], t=0..Pi, x=-sin(t)..sin(t), style=surface, color=yellow):
plots:-display(B1,B2,S1,S2, scaling=constrained, axes=normal, view=[-1.4..1.4,-1.4..1.4,-1.4..1.4]);

   
Addition. We can calculate the volume using a simple 1-integral if we notice that in the intersection with the plane  y=C  we have the square with the side  2*sqrt(1-C^2) :

V=2*int((2*sqrt(1-y^2))^2, y=0..1);

                                           V=16/3

Specify the value of A and remove one extra initial condition (should be 8).

See corrected file.

code_for_eigen_value_new.mw

Apply  ListTools:-Categorize  command to the list  convexhull(X), which will divide the entire list of points into equivalence classes according to the equality of the first coordinates. Then from each class choose the point with the smallest second coordinate.

A:=a^k*k/(a*a__0^k);
subs(a__0^(-k)=1/a__0^k, simplify(A));

Addition: here is a more general and programmatic way. Procedure  P  replaces any power with a negative exponent with a positive exponent.

P:=expr->evalindets(expr,`^`,t->`if`(sign(op(2,t))=-1,1/op(1,t)^(-op(2,t)),t)):

Examples of use:

P(simplify(a^k*k/(a*a__0^k)));  # The initial example
P(a^b-c^(-3)+1/(2^(-a)+c^(-b))+a^(-3)*b^(-3));  # More complicated example

The shorter ways:

eval([x,y,z], solve({2*x+y-3*z-3=0,x=-1+3*t,y=1+t,z=2-t}));
# or
solve({2*x+y-3*z-3=0,x=-1+3*t,y=1+t,z=2-t})[2..-1];
                                         [2, 2, 1]
                                {x = 2, y = 2, z = 1}

Another way is to use  geom3d :

with(geom3d):
intersection(P, line(l, [-1+3*t, 1+t, 2-t], t), plane(p, 2*x+y-3*z-3 = 0, [x, y, z])):
coordinates(P);
                                          [2, 2, 1]

Edit.

The curves  x(t)  and  y(t)  are periodic with some period T. If I correctly understood the question, then we should depict these curves on the cylinder of radius  T/(2*Pi) .

Example:

Eq:=diff(x(t),t)=4*x(t)-3*x(t)*y(t), diff(y(t),t)=x(t)*y(t)-2*y(t):
inc:=x(0)=1.5, y(0)=1:
sol:=dsolve({Eq, inc}, {x(t),y(t)}, numeric):
A:=plots:-odeplot(sol, [[t,x(t)], [t,y(t)]], t=0..5, color=[blue,red], thickness=3,  view=0..3.5):  
A;  # x(t) and y(t) in the range t=0..5

      
f:=s->eval(y(t),sol(s)):
T:=fsolve(f-1, 2..2.5);  # The period
r:=T/2/Pi;  # The radius of the cylindar
B:=plots:-odeplot(sol, [[t,x(t)], [t,y(t)]], t=0..T, color=[blue,red], thickness=3,  view=0..3.5):
B;  # x(t) and y(t) in the range t=0..T

F:=plottools:-transform((x,y)->[r*cos(x/r),r*sin(x/r),y]):
plots:-display(plot3d([r*cos(t),r*sin(t),h],t=0..2*Pi,h=0..3.5, style=surface, color=pink), F(B), axes=normal, view=[-0.43..0.43,-0.43..0.43,0..3.7]);  # The curves x(t) and y(t) on the cylindar

Curves_on_cylindar.mw

1. I do not know how to strictly prove that at  theta=Pi  there will be a singularity, but probably this is exactly so. If we set Digits:=100, we can calculate for  theta=Pi-10^(-95). We see that there is no stabilization, the values of the function  mu[1](theta)  are steadily decreasing.

2. I corrected your second file also.

1_new.mw

2_new.mw

f:= w^8+w^6+4*w^4+w^2+1: 
g:= w^16+2*w^14+9*w^12-2*w^10+44*w^8-2*w^6+9*w^4+2*w^2+1:
A:=[fsolve(f, complex)];
B:=[fsolve(g, complex)];
pA:=[Re,Im]~(A);
pB:=[Re,Im]~(B);
plots:-display(plot([pA, pB], style=point, color=[blue,red], symbol=solidcircle, symbolsize=15), plot([cos(t),sin(t),t=0..2*Pi], color=black), scaling=constrained); 
 

We get a better plotting quality if we use a different method. Also this method does not use  plots:-shadebetween  command, so it is suitable for older versions of Maple.
We use the fact that the projection of the solid on xOy plane is an ellipse.

restart;
eliminate({x^2+y^2-z^2=0, x+2*z=3}, z);
Student:-Precalculus:-CompleteSquare(%[2][],x);
((%=0)+12)/12; # The equation of the ellipse
a:=2: b:=sqrt(3): 
A:=plot3d(eval([x,y,sqrt(x^2+y^2)],[x=a*r*cos(t)-1,y=b*r*sin(t)]), r=0..1,t=0..2*Pi):
B:=plot3d(eval([x,y,(3-x)/2],[x=a*r*cos(t)-1,y=b*r*sin(t)]), r=0..1,t=0..2*Pi):
plots:-display(A,B, axes=normal, scaling=constrained);

p1 and p2 are the same plane, because the coefficients are proportional. So angle=0 degrees

I do not think that there is a simple formula expressing this integral for arbitrary numbers  a  and  m. It is reasonable to write this down as a procedure:

P := (a, m)->int((x-a)^m/x, x);

Examples of use:

P(2, 3); 
P(-2, 3); 
P(2, -3);

A:=a*diff(F(x,y),x,x)+b*diff(G(x,y),x,x)+c*diff(F(x,y),x,y)+d*diff(F(x,y),x,x)+e*diff(G(x,y),x,y)+f*diff(F(x,y),y,y)+g*diff(G(x,y),x,y);
sort(collect(A, diff), (x,y)->(has(x,F) and has(y,G)));

Edit.