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.

You're simply not correctly quoting the file name or including the file extension. Try this:

#Create list L and Matrix A before doing the save.

#Verify L is a list and A is a Matrix:
type([L,A], [list,Matrix]); #Should return true
save L, A, "mysavefile.mpl";
restart;
read "mysavefile.mpl": #End with a colon to avoid lengthy output.
type([L,A], [list,Matrix]);

If the last command returns true then it worked; if it returns false then report back here and we'll fix it.

I think that the above addresses all your questions, but to be explicit:

• Does the file need to be created first, before being saved to?

No, and if it has been created, the save command will simply overwrite it.

• Where should that file be placed? Desktop? Folder?

It makes no difference, as long as you save to and read from the same place. In the above example, I simply used whatever the default folder was.

• How do I read from the file after I have saved?

Use the read command.

• Can I save both L and A to the same file?

Yes.

• Do I need 2 commands?

No. If you want them in the same file, then you must save them with a single command.

I wrote a module to do it:

EnvelopePlot:= module()
export
   Extract:= (P::specfunc(PLOT))->
      [plottools:-getdata(P, "curve")][1][-1],
   
   Upper:= proc(C::Matrix)
   local X:= C[..,1], Y:= C[..,2];
      if not LinearAlgebra:-Equal(X, sort(X)) then
         error "CURVE isn't sorted wrt its 1st coordinate."
      end if;
      <X |
         CurveFitting:-ArrayInterpolation(
            C[select(
                 k-> Y[k-1] <= Y[k] and Y[k+1] <= Y[k],
                 [$2..op(2, [rtable_dims(Y)][1])-1]
              ), ..],
            X,
            method= spline
         )
      >
   end proc,

   Lower:= (C::Matrix)-> 
      < C[..,1] | -Upper(< C[..,1] | -C[..,2] >)[.., 2] >
;
end module:   

And here's an example of its use on an odeplot:

S:= dsolve({diff(f(x),x) = diff(cos(16*x)/(1+x^2), x), f(0)=1}, numeric);
P:= plots:-odeplot(
   S, [x,f(x)], x= -2*Pi..2*Pi, 
   numpoints= 2^10, color= black, legend= ["Signal"]
):
C:= EnvelopePlot:-Extract(P):
plots:-display(
   P, 
   plot(EnvelopePlot:-Upper(C), color= red, thickness= 3, legend= ["Upper envelope"]),
   plot(EnvelopePlot:-Lower(C), color= blue, thickness= 3, legend= ["Lower envelope"],
   axes= boxed)
);

The maximum number of objective-function evaluations is set by default to 50. You can change that by passing option evaluationlimit= 99 (or whatever number) to Maximize. I did this for your PD(10, 10, 10, 200), and there was no problem; the max number of evaluations that was actually used was 53.

Regarding the gap problem: I suspect that there's a bug in your code, but I'm not sure.

_U1 is a "bound", yet global, variable. It is the variable of integration. When Maple needs to make up a variable name to place in output, the name is of the form underscore-letter-integer. You can use subs to replace the name with another if you want. Note that t is used in the integral, so replacing _U1 with t would lead to an incorrect result. 

Z:= sqrt(1 - r^2*sin(theta)^2):
plots:-display(
   plot3d([r, theta, Z], r= 0..1, theta= -Pi..Pi, coords= cylindrical),
   plot3d([r, theta, -Z], r= 0..1, theta= -Pi..Pi, coords= cylindrical),
   plot3d(
      [1, theta, z], theta= -Pi..Pi, z= eval(-Z..Z, r= 1), coords= cylindrical, 
      style= wireframe, color= purple, thickness= 4
   ), 
   scaling= constrained, labels= [x,y,z]
);

#Volume integral in cylindrical coordinates:
Int(r, [z= -Z..Z, theta= -Pi..Pi, r= 0..1]):
% = value(%);

#Volume integral in cartesian coordinates:
Int(1, [z= -sqrt(1-y^2)..sqrt(1-y^2), y= -sqrt(1-x^2)..sqrt(1-x^2), x= -1..1]):
% = value(%);

It's clear by inspection that (x,y) = (0,0) is a solution for any t. It can be difficult for numerical methods to avoid converging on a trivial solution. Do you expect other solutions, and if so, can you give ranges in which you expect them to fall? They'll need to be ranges that don't include 0.

If we use a definite integral with a variable upper limit and some reasonable assumptions, then we'll get the obvious result based on the binomial expansion of (x-a)^m.

restart:
value(int((x-a)^m/x, x= a..x)) assuming x > a, a > 0, m::posint;
#Terms without x can be discarded as constants of integration:
select(has, %, x);
#Convert GAMMA to factorial:
convert(%, factorial);
subs(_k1= k, %);
simplify(%);

-x^m*(sum((-1)^k*a^k*x^(-k)/((-m+k)*k!*(m-k)!), k = 0 .. m-1))*m!+a^m*(-1)^m*ln(x)