Use  InertForm  package.

Example:

A:=<1,2; 3,4>:
B:=<3,4; 1,2>:
with(InertForm):
L := Parse("A*B"):
Display(L) = A.B;

A+`&lambda;I`=<1,2; 3,4>;  # Or

A+`&lambda; I`=<1,2; 3,4>;                     

To solve this problem,  combinat:-composition  command is useful. 

restart;
mondeg:=proc(var::list(name), d::nonnegint)
local n, L;
uses combinat, ListTools;
n:=nops(var);
L:=[seq(map(t->t-~1, Reverse(convert(composition(k+n, n),list)))[], k=0..d)];
map(t->`*`((var^~(t))[]), L);
end proc:

Example of use:

mondeg([x,y,z], 3);
         [1, x, y, z, x^2, x*y, x*z, y^2, y*z, z^2, x^3, x^2*y, x^2*z, x*y^2, x*y*z, x*z^2, y^3, y^2*z, y*z^2, z^3]

P:=6*x^2+x-2:
factor(P);
map(t->lcoeff(t)*``(t/lcoeff(t)), factor(P));
                                                             (3*x+2)*(2*x-1)
                                                            6*(x+2/3)*(x-1/2)

Does this match what you wanted?

For each n, I took 5 frames to lengthen the whole process:

restart;
with(plots):
f := x -> x;
L := 1: 
fs:=n->add(2*(-1)^(1+i)*sin(Pi*i*x)/Pi/i, i=1..n):
F:=n->plot([fs(n),f(x)], x=-L..L, color=[blue,red]):
display(seq(F(n)$5, n=1..25), insequence);

      


An alternative way to create an animation is to use  plots:-animate  command. Here is the code for your example:

restart;
with(plots):
f := x -> x;
L := 1: 
fs:=n->add(2*(-1)^(1+i)*sin(Pi*i*x)/Pi/i, i=1..n):
F:=n->plot([fs(n),f(x)], x=-L..L, color=[blue,red]):
animate(display,['F'(floor(n))], n=1..25, frames=100);

This method is preferable if the animation parameter changes in some continuous interval  a .. b .

Edit.

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