1. In the symbolic solution of a cubic equation, there is usually the imaginary unit  I, even if all the roots are real numbers (the so-called irreducible case). For example:

solve(x^3-3*x+1=0);
fsolve(x^3-3*x+1=0);

2. In general case, it is unlikely that something will be simplified. In special cases in addition simplify command you can also try  radnormal command.

3. In fact, this is the real number:
evalf(0.1159017153e-1*(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+3.632452480/(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+.6560037088+(1/2*I)*sqrt(3)*(-0.2318034307e-1*(820.0133-(3168.172220*I)*sqrt(3))^(1/3)+7.264904960/(820.0133-(3168.172220*I)*sqrt(3))^(1/3)));

                       0.6965691596-5.660254040*10^(-11)*I

The imaginary unit appears as a result of rounding errors.

In any case, you somehow have to refer to these rows. This is convenient if you use end-to-end indexing.

For your example:

row[1]:=[1,2,3];  
row[2]:=[4,5,6,7];
row[3]:=['red','blue',7];
Matrix(convert(row, list));

Note that Maple automatically aligns the lengths of the rows by adding zeros.

If you want to keep empty spaces in case of different lengths of lists, you can do this:
row[1]:=[1,2,3]:  
row[2]:=[4,5,6,7]:
row[3]:=['red','blue',7]:
m:=max(seq(nops(row[k]), k=1..3));
Matrix([seq([row[i][], ``$(m-nops(row[i]))], i=1..3)]);

Edit.

Here is another solution to the initial problem, based on the application of the Green formula. It requires parameterization of the boundary of the region. To do this, we first find the coordinates of the points  A, B, C, E, opening the modules in the corresponding regions. In the figure below, these regionss are highlighted with dashed green lines:

with(plots):
Eq1:=abs(abs(x-y)-abs(x+y))=2*y-x+1:
Eq2:=(x+1)^2+(y+1)^2=2:
BrokenLine:=implicitplot(Eq1, x=-2.5..1.5, y=-3..1.5, color=red, gridrefine=3):
Circle:=implicitplot(Eq2, x=-2.5..1.5, y=-3..1.5, color=blue, gridrefine=3):
Lines:=implicitplot([x-y=0, x+y=0], x=-2.5..1.5, y=-3..1.5, linestyle=3, color=green, gridrefine=3):
Points:=textplot([[0.18,-0.32,"A"],[-0.07,-0.55,"B"],[-0.4,-0.25,"C"],[-2.35,-0.95,"E"]], font=[roman,14]):
display(BrokenLine,Circle,Lines,Points, scaling=constrained);

By the code below we find the points A, B, C, E  and also the values of the parameter  t  in the parametric equation of the circle corresponding to the points  A  and  E :

solve({Eq1,Eq2}, explicit, useassumptions) assuming x+y<0, x-y>0, x>0;
A:=rhs~(convert(%, list));
solve({Eq1,x=0}, explicit);
B:=convert(rhs~(%), list);
EC:=simplify(Eq1) assuming x+y<0, x-y<0, x<0, y<0;
CB:=simplify(Eq1) assuming x+y<0, x-y>0, x<0, y<0;
C:=rhs~(convert(solve({EC,CB}), list));
solve({EC,Eq2}, explicit, useassumptions) assuming x+y<0, x-y<0, x<0, y<0;
E:=rhs~(convert(%, list));
Eq21:=[x=sqrt(2)*cos(t)-1,y=sqrt(2)*sin(t)-1];
t1:=rhs(solve(eval(Eq21, [x,y]=~E))[]);
t2:=rhs(solve(eval(Eq21, [x,y]=~A))[]);

                          A := [-7/13+(2/13)*sqrt(22), -17/13+(3/13)*sqrt(22)]
                                                     B := [0, -1/2]
                                                   C := [-1/3, -1/3]
                           E := [-19/17-(4/17)*sqrt(30), -9/17-(1/17)*sqrt(30)]

                   t1 := arctan((-(1/34)*sqrt(30)+4/17)/(-(2/17)*sqrt(30)-1/17))+Pi
                      t2 := arctan(((3/26)*sqrt(22)-2/13)/((1/13)*sqrt(22)+3/13))

Now everything is ready to get the final result. We use the procedure  Area  from here:

L:=[[A,B,C,E], [rhs~(Eq21), t=t1..t2+2*Pi]];
Result:=radnormal(Area(L));
evalf(Result);

   Result := (1/17)*sqrt(15)*sqrt(2)-(1/13)*sqrt(11)*sqrt(2)+Pi+arctan(2530/707-(479/707)*sqrt(11)*sqrt(2)+(1193/2121)*sqrt(15)*sqrt(2)-(440/2121)*sqrt(11)*sqrt(15))+419/663

                                              4.419820812

The result is identical to the result obtained earlier. Of course, the considered method is less elegant than the previous one, but it can be useful if  solve  command can not split the initial region into simple subregions. For example, it fails in this calculation (the intersection of three circles), although the expressions are polynomial:

rel:={(x-4)^2+y^2<=25, x^2+(y-3)^2>=9, (x+sqrt(7))^2+y^2>=16}:
solve(rel, [x,y]);
Warning, solutions may have been lost

The calculation of the area of region  rel can also be found here.

Area.mw
 


                          

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