X:= <21.2, 24.7, 20.8, 20.8, 20.3>:
Y:= <16.6, 19.7, 16.4, 16.8, 16.9>:
L:= Statistics:-LinearFit([1,x], <X|Y>, x);
plot([<X|Y>, L], x= min(X)..max(X), style= [point, line]);

The result returned by plottools[getdata] is a list, one of whose members is the matrix that you want. You need to extract the matrix from the list. Since the matrix is the third item in the list, that can be done by indexing:
M:= plottools[getdata](P1)[3];

Now M is what you referred to as a "simple" matrix, which you apparenty already know how to export.

IMO, there's only one good way to handle this very common situation. Make the first statement of the procedure:
if not x::realcons then return 'procname'(args) end if;
This is called returning unevaluated. All the attempts to correct this problem with quotes or is are crude and unstable workarounds.

Then, you should also apply evalf to x when it is being checked in an inequality.

What you want is as simple as seq(a__||k, k= 1..3), or, even better, a__||(1..3).

Here's a slick technique---one of a great many possibilities:

p:= randpoly(x):
select(type, allvalues~([solve(p)]), positive &under evalf);

Of course, your statement "All coefficients/variables in the polynomial are positive" is nonsense. 

I believe that the following handles equality testing for all root-of-unity cases, as shown by the subsequent verification:

restart:
RoU:= (k,n)-> exp(2*Pi*k*I/n):
RoU_Equal:= (a,b)-> evalb(evala(Norm(convert(a-b, RootOf))) = 0):

Verification code:

Test_Equal:= (a,b,c,n)-> 
   not (
      (RoU_Equal(RoU(a,n)*RoU(b,n), RoU(c,n)) xor irem(a+b-c, n) = 0) or
      (RoU_Equal(RoU(a,n)/RoU(b,n), RoU(c,n)) xor irem(a-b-c, n) = 0) or
      (RoU_Equal(RoU(a,n)^b, RoU(c,n)) xor irem(a*b-c, n) = 0) or
      (RoU_Equal(RoU(a,n)^(-b), RoU(c,n)) xor irem(-a*b-c, n) = 0)
   )
:

The following verifies every possible case, both true and false, with every possible base up to and including 9:

V:= <''a'', ''b'' , ''c'' , ''n''>:
for v in Iterator:-CartesianProduct([$0..8] $ 3, [$2..9]) do
    assign~(V=~v);
    if not Test_Equal(a,b,c,n) then
       error "Failure at %1, %2, %3, %4", a, b, c, n
    end if
end do:

Equal(<a>, <b>) is equivalent to evalb(a = b). I think that you have the mistaken idea that Equal is a more-sophisticated comparator. One that is actually more sophisticated is is(a = b).

Here's another recursive procedure that solves the problem.

bit:= proc(n::nonnegint, k::nonnegint)
local
   Put:= proc(T,n,v) T[n]:= v; T end proc,
   recur:= proc(n,k)
   local j; 
      if n=k then {table([seq(j=1, j= 1..n)])}
      elif k=0 then {table([seq(j=0, j= 1..n)])}
      else map(Put, thisproc(n-1, k), n, 0) union map(Put, thisproc(n-1, k-1), n, 1)
      end if
   end proc
;
   convert~(recur(n,k), list)
end proc:       

The time and memory values on the status line are only updated after a garbage collection. If the computation is at a stage where no garbage is being generated, then there'll be no garbage collection and no status updates.

If you're not familiar with the computer-science lingo, garbage is memory that is no longer accessible. For example, if you did

A:= [1,2];  A:= [3,4];

then the memory holding [1,2] would become garbage. Garbage collection is the process of collecting these scraps of unused memory and rearranging them into usable chunks of available memory. Garbage collection only occurs after a certain amount of garbage has piled up, typically many megabytes.

When you replace a list element with NULL, it changes the length of the list:  it is reduced by one. This behavior is different for Vectors: NULL can be a Vector element. 

I think that you can achieve what you want without a lot of intricate work with tildes. The procedure Vectorize, below, is essentially a form of unapply that takes a symbolic procedure and returns a procedure that's designed to work efficiently on float[8] rtables but will also work on other rtables or container objects. I believe that this completely automates the operation that you described.

Vectorize:= P-> proc(X) try map[evalhf](P,X) catch: P~(X) end try end proc:
A:= LinearAlgebra:-RandomMatrix((3000,3000), datatype= float[8]):
F:= Vectorize(D(x-> sin(x^2))):
FA:= CodeTools:-Usage(F(A));

memory used=0.54GiB, alloc change=68.67MiB, cpu time=4.33s, real time=3.67s, gc time=1.48s

The corresponding input is Eval(H1, t= 0.45). That's with a capital E so that it's inert. To activate the evaluation, you'd use value(%) (if you were actually running Maple, of course).

Use the -q (for quiet) flag on the cmaple command, for example

cmaple -q < myinputfile.mpl

In your module, include a procedure named ModuleLoad (you must use precisely this name). It may be declared local or export, but it needs to be one of these. Inside this procedure, include the statement :-Clr:= B. Right before the end module, include the statement ModuleLoad().

I don't understand why you want Clr to be global rather than a module export.

You need to use eval rather than subs. It is a lot like subs, but much more mathematically sophisticated, and it understands how to handle derivatives, integrals, etc. In constrat, subs is intended primarily for low-level programming. It just blindly substitutes what you tell it to.

Here's how to use eval in your situation:

g:= (x,y)-> diff(u1(x, y), x, x) + diff(u2(x, y), x, y):
f:= (y)-> eval(g(x,Y), [x= 0, Y= y]);

Note the distinction that I have made between lowercase y and uppercase Y.

Here's a completely different way:

g:= D[1,1](u1) + D[1,2](u2):
f:= (y)-> g(0,y):