Here's how to do it. I've coded it all as top-level code rather than as a procedure that produces each iteration. I removed all the underscores from your notation just to spare myself from typing them.

restart:

#Step 0:
B1:= A1+A2:  B2:= A1+A3:  B3:= A3+A4:  B4:= A2+A4:
C1:= `+`(A||(1..4)):
R:= [A||(1..4), B||(1..4), C1]:

#Step 1:
S01:= subs(seq(A||i= A||i||1, i= 1..4), R):
S02:= subs(seq(A||i= A||i||2, i= 1..4), R):
S1:= {S01[], S02[], zip(`+`, S01, S02)[]}:
nops(S1); #Confirmation
     27

#Step 2:
z:= P-> `+`(P[]) mod 2:
S2:= map(z, combinat:-choose(S1,2)) minus S1:
nops(S2); #Confirmation
     162

I adapted a cartesian-product procedure from the MaplePrimes archives:

CartProdSeq:= proc(L::seq(set))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    eval({subs(Seq=seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))});
end proc:

And with that machinery built, your main loop is simply 

Union:= S1:
for n from 3 while S||(n-1) <> {} do
     Union:= Union union S||(n-1);
     S||n:= map(z, CartProdSeq(S1,S||(n-1))) minus Union
end do;

I omitted the lengthy output. But S3 has 66 elements and S4 is empty.

surfdata requires that an Array of [x,y,z] data be 3D.  Since your M is a 441 x 3 Array, and since 441 = 21^2, I assumed that M corresponded to a 21 x 21 MESH.

M:= ExcelTools:-Import("C:/Users/Carl/desktop/Mdata.xlsx"):
M3:= ArrayTools:-Alias(M, [1..21, 1..21, 1..3]):
plots:-surfdata(M3);

If you evaluate r(X) at randomly chosen numeric values of X, then the value (i.e., the rank) will "almost always" be maximal. To quantify "almost alawys", use random integer X and compute the rank over a large finite field, where the probability that the rank of the numeric matrix is not maximal is a very small but easily computable number. Repeat independent trials until the product of the probabilities is below any epsilon that you want. Use enough different fields to lift the results to the reals.

Make f this:

f:= proc(t,M)
local n;
     if not [args]::list(numeric) then 'procname'(args)
     else evalhf(add(-4/Pi/(2*n+1)*sin((2*n+1)*t), n= 0..M))
     end if
end proc:

This involves passing M to f in addition to t. Then the plot call looks this:

plot(f(t,100), t= 0..2*Pi);

This runs in less than 1/10 second, using M=100.

Issue the command

Units:-UseSystem('CGS');

Names are by default symbols in Maple; you don't need to do anything to make them symbols. They stay symbols until you assign something to them. Your problem is caused by a lack of assumptions on your constants.

You do not need to enter with(Statistics) twice. You should not modify the global sequence constants. That sequence is for constants that have a definite value, like Pi.

restart:
with(Statistics):
h:= RandomVariable(Exponential(H)):
g:= RandomVariable(Exponential(G)):
assume(Pr > 0, Ps > 0):
PDF(Ps*g*h/(Pr*g+2), t);

This is the essence of what Acer was referring you to. I have vastly simplified it because I am not trying to emulate Mathematica or use this with procedures taking multiple arguments. We define a new binary infix operator which applies its right argument as a function to its left argument.

`&//`:= (A,B)-> B(A):

It can be used with % like this:

(x^2 - 1)/(x - 1); % &// simplify:

Or it can be applied directly to the left argument:

((x^2-1)/(x-1)) &// simplify;

But note carefully that in this case the extra parentheses on the left are required.

A chain of operations can be applied from the right without needing parentheses:

% &// f &// g;

Do you understand it now?

eq:= (4*y^2 + 8*y + 8*sin(y))/(y^2 +1) = 0:
numer(lhs(eq));

content(%);

4

%%/% = 0;

I don't know what criterion you used to select 8 as the coefficient to divide by: the largest coefficient? the most common coefficient? content finds the GCD of the coefficients. If you let me know your criterion, I can give you more code.

Include equation z=t and solve for [x,y,z]:

solve(
     [   x +  3*y +  5*z = 0,
       7*x +  9*y + 11*z = 0,
      13*x + 15*y + 17*z = 0,   
      z = t
     ],[x,y,z]
);

The command changevar is in the package student. You need to mention the package name in your command:

student[changevar](q, g, t);

Remove the with(linalg); that command is deprecated. In modern Maple, Jacobian is in the VectorCalculus package. Simply replace your Jacobian with VectorCalculus:-Jacobian. The result will be a 12x12 Matrix J whose entries can be accessed as J[1,1], etc.

No, you cannot access the member location inside a map invocation.

The following uses seq, but I'd be surprised if you could beat it timewise:

rtable([seq(f(i,M[i]), i= 1..1000)], subtype= rtable_options(M, subtype));

If f is thread safe (see ?threadsafe ), try using Threads:-Seq instead of seq.

Please post your code where seq is slowing down the computation. I find your statement about seq surprising. Perhaps someone here can improve the code.

There are two choices: Either you have Maple cancel units when and only when you ask it to, or you have it done automatically as expressions are entered.

To implement the former, use combine(..., units):

restart:
r:= 5*Unit(m/s):
t:= 10*Unit(s):
r*t;

combine(%, units);

To have it done automatically, you need to load a Units package environment. There are two environments that I am aware of: Natural and Standard. (Note that these terms, Natural and Standard, have nothing to do with the system of units being used. Rather, they refer to the method by which the units are entered into an expression.)

restart:
with(Units:-Natural):
r:= 5*(m/s);

t:= 10*s:
r*t;

restart:
with(Units:-Standard):
r:= 5*Unit(m/s):
t:= 10*Unit(s):
r*t;

Optimization:-Maximize can get stuck on a local maximum. The DirectSearch package (which must be downloaded from the Maple Applications Center ) has global optimizers. While Optimization works fine for this simplified example, it may not for your more-general case.

Sol:= dsolve(Sys, numeric, output= listprocedure):
y2:= eval(y2(t), Sol):  y3:= eval(y3(t), Sol):
# Above two lines same as Kitonum's
DirectSearch:-Search(abs(y2(t)-y3(t)), {t= 0..100}, maximize);

Use command StringTools:-WrapText.