This behavior is not specific to matrices. The object created by making assignments to a subscripted (a.k.a. indexed) name is always a "table". If a list (or Vector, Array, etc.) already exists, then the entries can be accessed by subscripting.  The most common way to create a list is to create a sequence and enclose it in square brackets ("[ ]"). For your example:

M:= [seq(Matrix(2,2, [[1, i], [i, i^2]]), i= 1..4) ]

Function application (f(...)) automatically (i.e. inherently to the syntax) distributes over list (or set ({})) creation, so the above could be

M:= [seq] (Matrix(2,2, [[1, i], [i, i^2]]), i= 1..4)

which I much prefer the look of, especially when the seq command spans several lines of code.

One may be tempted to create a list by repetively appending to an existing list, as in

M:= [];

for i to 4 do

    M:= [M[], Matrix(...)]

od;

This does create M as a list, but the problem is that it creates and discards 4 temporary copies of M.  The time complexity for this will be O(n²) where n is the list length. This will be a problem if n is large.  If one must use a for loop to create the list's entries, it is best to store them in a table within the loop, then convert the table to list outside the loop:

for i to 4 do  M[i]:= Matrix(...)  od;

M:= [seq](M[i], i= 1..4);

[1,2] + [3,4] *~ 3 *~ (x-1)^2;

Enter your integral in unevaluated form, i.e. as Int with a capital "I", and with symbolic e_0 as the upper limit:

t:= Int((1+121/304*e^2)^(1181/2299)*(1-e^2)^(-3/2)*e^(29/19), e= 0 .. e_0);

Then, simply,

plot(t, e_0= 0 .. 0.9);

Note to other readers: In this case, I needed to make sure that Maple did not try to evaluate the integral symbolically.  There is a very serious bug such that evaluating the integral (i.e. value(t)) causes the Maple kernel to utterly and completely crash after using about 11 seconds of processor time.

plots:-display(

   [seq](

      plots:-odeplot(

         dsolve(eval(sys, n= p), numeric)

        ,[x,g(x)]

        ,0..12

      )

     ,p in [0.1e-1, .1, .7, 1, 10]

   )

);

You could use the spreadsheet package.  For example,

> with(Spread):

> CreateSpreadsheet("mySS");

> M:= <<a,b>|<c,d>>:

> SetMatrix("mySS",M):

For example,

plot(1+x, x= -2..2, linestyle= dash, thickness= 3, color= green);

The order of the options does not matter.  For more details see ?plot,options ---probably my most frequently visited help page.

Apparently, fsolve has trouble with double integrals presented in the Int(..., [x= ..., y= ...]) form as opposed to the Int(Int(..., x= ...), y= ...) form.  Note the difference in the following two examples. 

restart:
J:= Int(a*x+y, [x= 0..1, y= 0..1]):
fsolve(J=0, a);
Error, (in fsolve) {x, y} are in the equation, and are not solved for

restart;
J:= Int(Int(a*x+y, x= 0..1), y= 0..1):
fsolve(J=0, a);
                          -1.000000000

Yet, for both forms type(J, freeof(x)) returns true. 

You can use the coeff command on either the series form or on the polynomial form that you had before converting it to a series.  If S is your series, then coeff(S,xi,0) and coeff(S,xi,2) will extract the expressions in your example.