One possible such command is codegen[cost], but in this case the exponentiations (with small integer exponents) would be counted as an appropriate number of multiplications, because that's how the expression would be numerically evaluated.

There's no need to use arrow expressions; they just complicate things in this case. The following code does what you want. It works for any number of input expressions. It doesn't rely on any specific names being used such as x, y, or t.

linearise:= proc(S::list(algebraic), P::list(anyfunc(name)=algebraic))
local F:= lhs~(P), t:= op~({F[]})[];
   if nops([t]) > 1 then 
      error "Expected exactly 1 independent variable but received", t
   end if;
   diff~(F,t) =~ thaw(map(mtaylor, subs(F=~ freeze~(F), [S,P])[], 2))
end proc:

#Example usage
a:= x(t)*(y(t)-1);
b:= 4 - y(t)^2 - x(t)^2;
linearise([a,b], [x(t)= 0, y(t)= -2]);

Here's another procedure. This'll work in any number of dimensions, as long as the hyperfaces are parallel to the coordinate hyperplanes.

Vertices:= proc(P::list, Q::list) 
local x,p;
   seq([seq(x, x= p)], p= Iterator:-CartesianProduct(zip(`[]`,P,Q)[]))
end proc:

Vertices([1,1,2], [3,4,5]);

I don't know whether it can be done with geom3d.

I suggest that you stop using ancient and deprecated stats package and instead use the newer Statistics package. Then

x:= Statistics:-Sample(Uniform(0,1), 10);
add(x)^a;

Since the above is hardly worth doing, I think that what you actually want is

add(x_i^a, x_i= x);

Right-click on the Matrix or Vector output and select Browse.

The following does exactly what you're trying to do using a different ODE (because you didn't supply yours) and using a list of three operators as the first argument to plot3d. Because of the option remember, td is only called once (instead of three times) for each mesh point of the plot.

ans:= dsolve({diff(y(t),t$2)=y(t), y(0)=a, D(y)(0)=0}, numeric, parameters= [a]):
td:= proc(x,tt)
option remember;
   ans(parameters= [a= x]);
   eval([t, y(t), diff(y(t),t)], ans(tt))
end proc:

TD:= (k::{1,2,3})-> (x,tt)-> td(x,tt)[k]:
plot3d(TD~([1,2,3]), 0..3, 0..3, labels= [t, y, `y'`]);

In Maple, several binary infix operators {=, <, .., ^, ::, etc.} always have exactly two operands. For =, each of the operands can be anything---including NULL, a sequence, or another equation---as long as it's enclosed in parentheses. So (A = B) = C and A = (B = C) are both allowed. I realize that the display of the parentheses in printed output is not ideal.

This is implemented in the kernel, so it works regardless of which GUI (Standard, Classic, command-line) or GUI mode (document, extended typesetting, etc.) you use.

P:= proc(u::{function, indexed}, shift::anything:= 1)
local pos:= `if`(procname::indexed(posint), op(procname), 1);
   subsop(pos= op(pos, u) + shift, u)
end proc:

P[1](u(i,j));
                          u(i + 1, j)
P[2](u[i,j,k], 3);
                          u[i, j + 3, k]

This allows you to have a single operator that does it all. It doesn't care whether you've used square brackets or parentheses. It doesn't care how many indices there are. And you can change the shift amount from the default of 1,

This can be prevented by using forget(fsolve) or gc() between the fsolves. An automatically occuring garbage collection would also clear the cache. Anytime that you see the status bar update, there's been a garbage collection.

Another approach is to simply use fsolve(f) in both cases. This'll work with or without quotes.

Yes, the problem is likely high memory usage due to large symbolic expressions. We call it expression swell. A direct symbolic solution is not easily amenable to distributed (multi-core) processing. Nor is any memory-bound process a good candidate for multi-core: Distributed processing increases memory usage significantly (which isn't a problem if the process is distributed over a network of computers).

Depending on the complexity of your symbolic coefficients, it may be possible to interpolate an exact symbolic solution from exact numeric solutions of the system numerically evaluated at numerous values of the parameters. This technique would be amenable to distributed processing. But first think about whether it's worth doing vis-a-vis my Reply above.

plot([seq([i, eval(u[i,20], Sol[])], i= 1..25)], style= point);

You can change the symbol used for the points, their size (symbolsize), and their color:

plot(
   [seq([i, eval(u[i,20], Sol[])], i= 1..25)],
   style= point, symbol= solidcircle, symbolsize= 16, color= red
);

When your ODEs involve a function that can only be given by a numeric procedure, then you need to use option known with dsolve, like this:

f:= proc(x,y)
   if [x,y]::list(numeric) then 
      x^2+y^2; #My numeric code
      #Put your numeric code here.
   else
      'procname'(args)
   end if
end proc:
 
Sol:= dsolve({diff(y(x),x) = f(x,y(x)), y(0) = 1}, numeric, known= [f]);

The plot is not empty. Looking at the plot structure shows that it's full of white arrows. This is further confirmed by including the option arrows= THICK. What is shading= none supposed to do, in your opinion? It seems to me that white arrows is reasonable interpretation of "no shading".

You asked for comments about your code: Your "union" algorithm is fine---it certainly seems like the natural approach---but your Maple implementation is very slow due to a nuance of Maple: It is very inefficient to create large sets by adding elements one at a time because every modification of a set requires the entire set to be moved in memory. The soluton to this problem is usually to accumulate stuff temporarily in Maple tables and create the big set (or list) later. Here's code that does that.

mondeg2:= proc(V::list(name), d::nonnegint)
local k, v, mon:= table([0= {1}]);
   for k to d do mon[k]:= `union`(seq(v*~mon[k-1], v= V)) end do;
   op~([entries(mon, 'nolist')])
end proc: 

Notice also that it's much shorter.

The code using coeffs and expand that I posted earlier is faster than the above. I believe that that's because expand has been highly optimized.

The following worksheet has timings for all three.
 

Download mondeg.mw

restart:
Digits:= 20:
CodeTools:-Usage(plot(x-> exp(1e9*ln(1-x)), 0..1e-15));
memory used=7.76MiB, alloc change=32.00MiB, cpu time=63.00ms, real time=69.00ms, gc time=0ns