As for Kitonum's solution, the following procedure returns all occurences of the maximum value:

findMax := proc(M::Matrix)
   local m := max(convert(M,set));
   select(x -> rhs(x) = m,op(convert(M,table)));
end proc:

Using polynomial division:

c := rand(1..9):
l := [seq(c()*x + c(),i = 1..4)];
not `or`(seq(seq(Algebraic:-Divide(l[i],l[j]),j = i+1..4),i = 1..4));

PS: Note that this method discards not only lists containing factors differing by an overall sign (your item 3), but generally lists containing factors proportional to each other.

Using lists of lists as intermediates:

removeZeroRowsAndCols := (M::Matrix) ->
   convert(remove(hastype,convert(
   convert(remove(hastype,convert(
      M,
   listlist),list(0)),Matrix)^%T,
   listlist),list(0)),Matrix)^%T:

Concerning the error: If your computer runs Windows, then try using forward slash, /, as separator in your paths.

Elaborating on Mac Dude's suggestion on making a package: If your procedures can be naturally divided into some categories, then you might want to consider making a package for each such category. These (sub)packages can then be used in some main module (package), say, either in their long form, package:-procedure, or by loading them with the command uses, see the help page ?uses. This might provide you with a nice modularization of the overall code.

I am not sure whether this proposal of mine satisfies your needs, but perhaps you can use selectremove. An example:

expr := 2*x^2 + 2*y^2 + a;
s,r := selectremove(has,expr,{x,y});
convert(s,elsymfun) + r;

I think that setting scaling = constrained in pointplot would do the trick, see the help page ?plot options, in combination with some appropiate scaling of the y coordinate, for instance, 1/2*(x + 1 - i).

The following is an example of a module inside a module:

m1 := module()
   export m2;
   m2 := module()
      export n;
      n := (x) -> 2*x;
   end module;
end module:
m1:-m2:-n(10);

Is that what you are looking for?

At the classical level, using Grassmann numbers, you could define a Dirac spinor in Minkowski space, say, as follows:

with(Physics):
Setup(
   mathematicalnotation = true,
   coordinatesystems = (X = (t,x,y,z)),
   anticommutativeprefix = psi   # The Grassmann-odd property of the components
):
dirac := Vector(4,(a) -> psi||a(X));

Perhaps the randomize command will be useful:

with(RandomTools):
randomize():
w:=Generate(list(posint(range = 3),15));

Concerning question 1:

(m - ap1)^2;

Concerning question 2, and expression ap2:

ap2 := numer(ap2) / denom(ap2);

I do not understand why you want to collect with respect to m - T; it seems to me that there is not really much to collect in ap2. But perhaps I misunderstand what you want to achieve.

Concerning question 3, and expressions ap3 and ap4: I think ap3 should be written as

ap3 := diff(t1(T),T) = (m-T) / (m - t1(T));

with t1 inside diff being explicitly dependent on T. I have also changed Diff to diff, but that might not be in line with what you intend. Note that Diff is not the same as diff.

I am having trouble deciphering what you intend in ap4. Could you explain?

Concerning question 4: The csgn function is, loosely speaking, the sign function. You can use the assume facilicities of Maple to treat the positive and negative cases, respectively:

simplify(ap5) assuming m - t1 > 0;
simplify(ap5) assuming m - t1 < 0;

I can see that you use Diff in ap5, rather than diff. Is that intentional?

Concerning question 5: In order for Maple to simplify ap6, you have to let it know something about beta. If it is real-valued, you could do as follows:

ap6 := simplify(ap6) assuming beta::real;

which brings you some of the way to your by-hand derived expression.

Perhaps a lowtech solution, but time() might be useful. A trivial example:

t := time():
for i from 1 to 100000000 do
   if time() - t > 1. then error i end if
end do:

With a for loop:

eqs := {}:
for i from 1 to 3 do
   eqs := eqs union convert((A(W[i]) . Phi[i])[1..2],set);
end do:
eval(<Phi[1]|Phi[2]|Phi[3]>,solve(eqs));

Or without a for loop:

eval(
   <Phi[1]|Phi[2]|Phi[3]>,
   solve(`union`(seq(convert((A(W[i]) . Phi[i])[1..2],set),i = 1..3)))
);

The syntax is along the following lines:

m;   # mass
r;   # radius
J := m*r^2;
m[a];   # mass of part a
m[b];   # mass of part b

Note that comments must be placed after #'s.

Warning: See update below.

The problem is in reality a Lagrangian multiplier problem, see ?LagrangeMultipliers. But since I cannot get that working, I proceed in a somewhat less technical manner:

Step 1: Insert z = a + I*b into the condition, and perform evalc, thus (correctly) assuming that a and b are real numbers:

eq := abs(z + 3 - 2*I) + abs(z - 3 - 8*I) - 6*sqrt(2);
eq := evalc(eval(eq,z = a + I*b));

Step 2: Solve the equation to satisfy the condition:

sol := solve(eq,{a,b});

It appears that a = b - 5, b arbitrary.

Step 3: Insert that constraint into the modulus of z = a + I*b, and look for extrema by the standard method of differentiation:

m := evalc(eval(abs(a + I*b),sol));
solve(diff(m,b),{b});

There seems to be only a single extremum at b = 5/2 (I wonder, is that actually correct?), at which the modulus of z is 5/sqrt(2), this value being a local minimum because the second derivative of m at b = 5/sqrt(2) is positive.

Update: There is something wrong in step 2 above. For if the solution is reinserted into the condition, zero does not result:

f := simplify(eval(eq,sol));
plot(f,b);

Now I am confused!

Using the package VariationalCalculus, you can do as follows:

with(VariationalCalculus):
EulerLagrange(L,t,[z(t),v(t)]);

PS: I am somewhat confused about the h-indexed z(t) of yours, i.e., by z(t)[h]. Is that really intended?