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

D is reserved for derivatives, so I just changed your D to D1. I could theoretically use D as a variable, but I'd just rather not (ask for trouble).

restart:
E[1]:= A+B+C+D1:
E[2]:= A-B+C-D1:
E[3]:= C:
V:= indets({entries(E)}, name): #So, V = {A, B, C, D1}.
EQ:= x*E[1] + y*E[2] + z*E[3] = w*A:
#We form a system of 4 (nops(V)) equations by making 4 sets of random integer
#assignments to the entries of V, using those in EQ, asking for
#integer solutions with isolve, and setting the integer parameter Z to 1. 
S:= eval(isolve({'eval(EQ, V=~['rand()'$nops(V)])'$nops(V)}, {Z}), Z= 1);

The result returned by isolve is

     {w = 2, x = 1, y = 1, z = -2}

I think that that guarantees (setwise, not pairwise) relatively prime solutions, if they exist at all. If I'm wrong about that, it's trivial to divide out the gcd.

Note that the expansion of (x+y+z)^k contains every monomial of degree k, with some integer coefficients that we can just ignore. That idea leads to this very simple procedure:

mondeg:= proc(V::list(name), d::posint)
local T, k, P:= add(V);
   coeffs(expand(add(P^k, k= 0..d)), V, T);
   {T}
end proc:

The usage is just like yours and Kitonum's:

mondeg([x,y,z], 3);

     {1, x, y, z, x^2, x^3, y^2, y^3, z^2, z^3, x*y, x*y^2, x*z, x*z^2, x^2*y, x^2*z, y*z, y*z^2, y^2*z, x*y*z}

Update, to remove add and to simplify based on vv's +1 idea:

mondeg:= proc(V::list(name), d::nonnegint)
local T;
   coeffs(expand(`+`(V[],1)^d), V, T);
   {T}
end proc:

If you're using Maple's "Standard" GUI, then you can use Explore rather than display to produce animations that play immediately:

Explore(plot([f(x),fs(m)], x= -L..L), m= 1..5, autorun, loop);

You need to start the animation ("play" it) by using the context menu or by selecting it and using the controls on the lowest toolbar.

Instead of all that typing, why not simply enter
< 1, 2, 3, 0;
   4, 5, 6, 7;
   'red', 'blue', 7, 0
>;

?

Perhaps I don't understand your purpose. Are you trying to facilitate the GUI entry Matrices, or do you have some programmatic reason for wanting to do it the way that you suggest?

I'm not at all opposed to Maple's object modules; indeed, I was going to post a solution using them had not Joe beat me to it. But there's an easier way to implement multiple independent instantiation of objects (in the OOP sense): You simply need a procedure that returns the module. In the following, I've simply wrapped my module from my previous Answer with proc() ... end proc. That's literally all there is to it.

CountProc:= proc() 
   module()
   local 
      m:= 0,
      ModuleApply:= proc() m:= m+1 end proc
   ;
   end module
end proc:

Now you can create independent counters, as in

counter1:= CountProc(): counter2:= CountProc():
counter1(), counter2(), counter1(), counter1(), counter2();

     1, 1, 2, 3, 2

Some writers call such a procedure a "module factory".

If you want to overload operators, such as what Joe did with !, then you pretty much need to use option object. There's a rather cruddy way to do it without option object, but it's so cruddy that I don't think that it's worth mentioning.

If the procedure CountProc had parameters, then the module could refer to them. This can be a very powerful technique.

This is precisely what modules are for. In this case:

CountProc:= module()
local 
   m:= 0,
   ModuleApply:= proc() m:= m+1 end proc
;
end module:

Now,

CountProc(); CountProc(); CountProc();

returns exactly the same results that yours did, but m is local.

Essentially, a module is just a procedure with persistent local variables. Some of those local variables are usually themselves procedures. Some of the local variables may be designated as exports (they are declared with keyword export instead of keyword local). This means that they are directly visible and changeable from outside the module.

In my example, ModuleApply is a keyword that says that the module's name, CountProc, can be used as if it were a procedure; and, in that case, the procedure ModuleApply is what is actually called. You must use the keyword ModuleApply to do this.

If R is a list instead of a Vector (which is the way you presented it), do

(vals,Numbs):= (rhs~,lhs~)(select((x-> x::numeric and x > -2 and x < 3)@rhs, [op(2, <R>)[]]));

     vals, Numbs := [-1, 2, 1, -1, 2.2, .124], [2, 6, 8, 9, 10, 14]

If R is in fact already a Vector, just change <R> to R in the above.