It is doubtful that increasing the limit will provide a fix. If you post the model here I might be able to make a recommendation.

You can write a procedure that does what you want.  For example:

NonCommutativeMatrixProduct := proc(A :: Matrix, B :: Matrix, prod)
local i,j,k,m,n,p,q,M;

    (m,n) := op(1,A);
    (q,p) := op(1,B);
    if n <> q then
        error "incompatible matrix dimensions";
    end if;

    M := Matrix(m,p);

    for i to m do
        for j to p do
            M[i,j] := add(prod(A[i,k],B[k,j],_rest), k = 1..n);
        end do;
    end do;
    return M;
end proc:

A := Matrix(2,3,'symbol'=a):
B := Matrix(3,4,'symbol'=b):

NonCommutativeMatrixProduct(A,B,`&^`);

Using your example

with(Physics);
Setup(anticommutativeprefix = {Yc, y});
C[1] := Matrix([[X[1], Y[1]], [y[1], x[1]]]);
C[2] := Matrix([[X[2], Y[2]], [y[2], x[2]]]);
Physics[`.`](C[1], C[2]);
Physics[`.`](C[2], C[1]);

Physics[`.`](C[2], C[1]);

NonCommutativeMatrixProduct(C[1],C[2],Physics[`.`]);

Physics[`.`](C[2], C[1]);
               [X[1] X[2] + y[1] Y[2]    X[2] Y[1] + Y[2] x[1]]
               [                                              ]
               [y[2] X[1] + x[2] y[1]    Y[1] y[2] + x[1] x[2]]


NonCommutativeMatrixProduct(C[1],C[2],Physics[`.`]);
               [X[1] X[2] + Y[1] y[2]    X[1] Y[2] + Y[1] x[2]]
               [                                              ]
               [y[1] X[2] + x[1] y[2]    y[1] Y[2] + x[1] x[2]]

For newer releases of Maple you can assign a default value that is different from the declared type and use that to determine whether an argument has been passed. Here NULL might be appropriate

proc( ...  S :: symbol := NULL, L :: list := [], ... )
    if S = NULL then <do something> end if;

(**) y := sum(a[n]*x^(n-2)*(n-1)*n, n = 0 .. infinity);
                               infinity
                                -----
                                 \            (n - 2)
                          y :=    )     a[n] x        (n - 1) n
                                 /
                                -----
                                n = 0

(**) op(1,y);                                          
                                       (n - 2)
                                 a[n] x        (n - 1) n

The problem lies with Maple's 2D parser and that input.  Mucking around with your input I can eventually get Maple's parser to accept it and work, however, it isn't clear what is being changed.  For this and similar reasons I never use 2D input, instead I use Maple input, which avoids these weird issues.

Because only the nonzero entries are stored in the Matrix, you could use the following method to extract the set of indices of the nonzero elements (here 1's):

 map([lhs],op(2,M));

I doubt that Maple can compute the closed-form expression by directly manipulating the sums.  It is possible to use Maple to compute the result by explicitly evaluating the expressions to generate a sequence, using gfun to convert to a recurrence equation, and then using rsolve.

First define a procedure that evaluates w

w := proc(i)
    if i :: numeric then
        if i=0 then 1 else alpha end if;
    else
        'procname'(i);
    end if;
end proc:

Next define a procedure that evaluates the summations

y := proc(h)
option remember;
    add(add(add(w(i)*w(i+j), i = 0 .. l-1), j = 1 .. h-l), l = 1 .. h);
end proc:

A little experimentation reveals that y generally evaluates to an expression y1(h)*alpha + y2(h)*alpha. So

y1 := h -> coeff(y(h),alpha,1):
y2 := h -> coeff(y(h),alpha,2):

We can now solve for an explicit expression for y1:

gfun[listtorec]([seq(y1(h), h=0..30)],u(h));
Y1 := rsolve(%[1],u(h));
                                     Y1 := -2 h - 1 + (h + 1) (h/2 + 1)

Alas, the same method, unassisted, doesn't work for y2.  However, let's look at the generating function

data := [seq(y2(h), h=0..20)]:
gfun[guessgf](data,h):
factor(%[1]);
                                       3
                                      h
                                   --------
                                          4
                                   (h - 1)

The h^3 in the numerator acts as a shift.  So let's try shifting the data

gfun[listtorec](data[4..], u(h)): 
Y2d := rsolve(%[1], u(h));
                                (3 + h) (h + 2) (h + 1)
                        Y2d := -----------------------
                                          6

From that we can manually create a piecewise expression for y2

 Y2 := piecewise(h<3, 0, eval(Y2d,h=h-3));
                        {         0                  h < 3
                        {
                  Y2 := { h (h - 1) (h - 2)
                        { -----------------        otherwise
                        {         6

If we limit the range of h to nonnegative integers, then this is equivalent to

Y2 := eval(Y2d, h=h-3);
                                  h (h - 1) (h - 2)
                            Y2 := -----------------
                                          6

w := proc(i)
    if i :: numeric then
        if i=0 then 1 else alpha end if;
    else
        'procname'(i);
    end if;
end proc;

y := Sum(Sum(Sum(w(i)*w(i+j), i = 0 .. l-1), j = 1 .. h-l), l = 1 .. h):
value(eval(y, h=3));
                                             2
                               3 alpha + alpha

I recovered what was readily recoverable, but from the look of it that was just the write-up from the instructor.  That is, there was nothing recoverable past the following caveats:

  Remember to shrink your graphs and save your Maple file frequently as you go!
  Also remember to highlight or italicize your answers when typing in a text box.

Note that Maple periodically saves backups. See ?backup. Possibly you can recover a useful version from that?

You can work around the first limitation by replacing the output of Step(), which returns the Heaviside function, with an explicit piecewise:

ResponsePlot(sys, piecewise(t<0,0,1));

Alternatively you can do the equivalent (see ?Heaviside for details):

_EnvUseHeavisideAsUnitStep := true:
ResponsePlot(sys, convert(Step(),piecewise));

What version of Maple are you using? I'm not seeing the second error (maxfun exceeded). The way to work around that is to do

ResponsePlot(sys, Sine(), duration=40, dsolveargs=[maxfun=0]);

Assuming it exists, Maple can solve for the limit:

req := a(k) = (2/5)*a(k-1)*ln((1/10)*(exp(1)*x/a(k-1))^10)/ln(10^a(k-1)*(exp(1))^4):
# Assuming limit exists, then at limit, a(k) = a(k-1)
leq := eval(req, a=(x->a)):
simplify(leq,'symbolic'):
a_lim := solve(%,a);
                                           x
          a_lim := -----------------------------------------------------------
                    exp(LambertW(1/4 x exp(-1/10 ln(10)) ln(10)) + 1/10 ln(10))
evalf(eval(a_lim, x=1));
                                             0.5716077045

In case it isn't clear, you can also plot expressions of the state-variables.  That allows, say, plotting y vs x:

odeplot(p, [[x(t),y(t)]],-4..4);

Could you upload the .msim fle that has the error?

As specified, the differential equation has two branches.  The easiest thing to do is to rewrite it, forcing one branch.  There is also a problem with the initial conditions, you don't want to specify D(theta)(0), that is not allowed for a first order system.  Is the variable l supposed to be the parameter longueur?  Here is a rewrite that works, but you probably need to tweak it

tipe6 := proc (longueur, thetaInitial, rayonBoule, nombrePoints)
local eq, sol, valeurstheta, x, y, plotseq;
global g,l,t,theta;
    #eq := ((1/2)*longueur^2+(1/5)*rayonBoule^2)*(diff(theta(t), t))^2+g*l*(cos(theta(t))-cos(thetaInitial)) = 0;
    eq := diff(theta(t),t) = (-(10*g*l*cos(theta(t))-10*g*l*cos(thetaInitial))/(5*longueur^2+2*rayonBoule^2))^(1/2);
    sol := dsolve({theta(0) = thetaInitial, eq}, theta(t), numeric);
    valeurstheta := seq(rhs(sol((1/30)*i)[2]), i = 1 .. nombrePoints);
    x := seq(longueur*sin(valeurstheta[i]), i = 1 .. nombrePoints);
    y := seq(longueur*(1-cos(valeurstheta[i])), i = 1 .. nombrePoints);
    plotseq := [seq(plot([[x[i], y[i]]], 'style' = point), i = 1 .. nombrePoints)];
    plots:-display(plotseq, 'insequence' = true)
end proc:

g := 10: l := 1:
tipe6(1,Pi/4,10,100);

While the help page for the WheelAxle component states that the equation for torque is tau = r*f_T*dir, in fact it actually is tau = -r*f_T*dir. The choice of sign depends on how the wheel is defined.