The problem lies in the way you are using indexed variables.  Consider the simple example

g[j] := j^2:
p := j -> g[j]:
p(3);
                   g[3]

That returns g[3], and not the 9 that your scheme would require.  One way to fix this is to use functions instead of indexed names:

g := j -> j^2:
p := j -> g(j):
p(3);
                             9

Here is a simple workaround.  Use the postscript plot device and create the file.  Then manually edit it so that the S commands are changed to P commands (which fill the polygon).  You might have to play around with this depending on what else is plotted.  Here is what I did (this assumes you have access to the sed command on your system):

psfile := "plot.ps":
plotsetup(ps, plotoutput=psfile);
plot(sin, style=point, 'symbol=solidbox');
# Use sed to convert all S commands to P commands; save the original
# as a backup.
ssystem(sprintf("sed -i.bak 's/^S/P/' %s", psfile));

That generates solid squares.

If you don't have access to sed, you could do the conversion with Maple:

tmp := FileTools:-Text:-OpenTemporaryFile();
do 
    line := readline(psfile);
    if line = 0 then break; end if;
    line := StringTools:-RegSubs("^S$"="P",line);
    fprintf(tmp, "%s\n", line);
end do:
close(tmp);
FileTools:-Rename(tmp,psfile,'force');

The reason that you cannot use _dj01ajc is that internally fdiff is setting Digits to 17, and that integration method uses hardward floating points, so will not work if Digits is greater than evalhf(Digits).  One way to get a result is to reduce Digits to 8. 

Digits := 8:
fdiff(Tproc, [1], [1,2,3]);
fdiff(Tproc, [2], [2, 10, 2.5]);

The fdiff routine should be generating an error message rather than returning the result unevaluated; I'lll submit an SCR.

Hard to say without inspecting your code.  Can you upload it?

The problem lies in the assignment to F₁, you omitted a multiplication immediately following x⁵. Consequently Maple interprets it as a function, sort of.  That is the 2D parser makes x the function and raises that to the power of 5.

It would be helpful if you included the source of what you actually did.  The description is not precise.

roll := rand(0 .. 1):
s0 := 0: s1 := 0:
to 100 do
   s0 := s0 + roll();
   s1 := s1 + roll();
end do:

s0,s1;
                              50, 46

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