The command is

eliminate(eqns, {X,Y});

The solution may be messier than you were hoping for, with a few nasty RootOfs.

To run a .mpl file in your current Maple session, the command is simply

read "filename":

Your attempt using cmaple.exe would not run in the same Maple session (or kernel), so it would be of no use to you. However, it can occasionally be useful to run an external Maple session like you were trying. This is done with the system or ssystem commands.

A color function can be arbitrarily complicated on the inside, but its signature must be relatively simple: It must take two real arguments and return a single real result. It cannot return a list. If you want to use a 3D colorspace like HSV or RGB, then each coordinate color needs a separate function, like this:

restart:
(X,Y,Z):= (((u,v)-> sech(u)*cos(v)), ((u,v)-> sech(u)*sin(v)), ((u,v)-> u - tanh(u))):
F:= z-> z:
Ans:= proc(u,v) option remember; evalf(Re(F(v+I*exp(u)))) end proc:
`2*Pi`:= evalf(2*Pi):
ColorFunc||(H,S,V):= (
   ((u,v)-> `if`(Ans(u,v) > `2*Pi`, 0, Ans(u,v)/`2*Pi`)),
   ((u,v)-> `if`(Ans(u,v) > `2*Pi`, 0, 1)) $ 2
):
plot3d(
   [X,Y,Z], 0..3, 0..`2*Pi`, grid= [2^5 $ 2], 
   color= [ColorFunc||(H,S,V), colortype= HSV],
   lightmodel= none
);

In your specification of equ6, you've used F[1](0). This is logical, and Maple should be able to figure it out because you've specified the value of F [1](0) in Bcs21. Unfortunately, dsolve wasn't designed to take this into account. You can replace F [1](0) with a; or, a more-sophisticated treatment is 

S21:= dsolve({eval(equ6, {Bcs21}), Bcs21}, numeric);

There are many reasons in Maple why it's better to set the value of a parameter with eval rather than := (mostly, it improves the clarity of the exposition). You've stumbled upon one of them. This is doubly true when the parameter is inside an rtable (a Vector, Matrix, Array, or anything constructed with <...>). If you insist on using :=, there is a command to force the assignments to take effect inside of the rtables: rtable_eval. Thus,

z:= A*<x,y>:
A:= 1:
rtable_eval(z, inplace):
plots:-fieldplot(z, x= -1..1, y= -1..1);

Note carefully however that the rtable_eval(..., inplace) changes z from an rtable with a parameter to a purely numeric rtable. To retain the original z, don't use inplace:

restart:
z:= A*<x,y>:
A:= 1:
plots:-fieldplot(rtable_eval(z), x= -1..1, y= -1..1);

I think that it's unfortunate that Maple's rtables were designed to be used for two conflicting purposes: as containers for efficient computation, and as realizations of the mathematical concepts of vectors and matrices. The command rtable_eval exists as a bridge between those purposes. In the long run, I think that it's better to use lists to represent mathematical vectors, and listlists to represent mathematical matrices.

You can view a 3D plot as if it were a 2D plot by including the option orientation= [270,0].

Since n is a bound variable to each sum, the n in each sum is logically a different variable, although it's the same "name" (or address) as far as Maple is concerned. So, your question doesn't quite make sense to me. If the result could be expressed as a single summation, that'd be a different matter. Then you could say

eval(subs((n= 0..infinity)= (n= 2..2), %));

and get a meaningful answer. You can also give this command in the case that you presented---the syntax allows it---but as far as I'm concerned, the result is nonsense.

A meaningful question that you could ask is What is the xi^2 term in the resulting expression? This question can be answered by

coeff(eval(%, infinity= 2), xi, 2)*xi^2;

You need two sets of unevaluation quotes around each name in KeepThese:

one:= 1: two:= 2: three:= 3:
KeepThese:= {''one'',''two''}:
unassign~({anames('user')} minus KeepThese);

Why not simply find the zeros of the Bessel function and then divide by R?

Whatever algorithm you're using does not allow complex numbers, yet they are being generated. Did you perhaps try to use capital I as a variable? That's not allowed; capital I is the complex unit sqrt(-1).

If that doesn't help, then you'll need to post what you entered to get to the above error message.

You asked:

• It seems me that this technique enable to adapt the procedure  (DO_CONSTRAINT) to different types of inputs that we can try on it (lists or scalars). Do I right ?

Yes.

• May you give me more informations on how it works ?

Since the code seems self-evident to me, you'll need to ask a more-specific question about what you don't understand. If you're wondering whether there's anything special about the name _DO_CONSTRAINT, there isn't; the procedure could just as well be named foo.

• Is it the only possibility to do this ? In other words, is there a possibility to make a similar code and also concise but without using this local variable ?

It is a common operation to adapt a procedure to accept different forms of input. There are many ways to do it. The way that you show seems crude in my opinion. I think that the most common way to do this is to define the procedure for basic input and then call it recursively for other input. So only one procedure is required. Like this:

DO_CONSTRAINT:= proc(phi::{list, scalar}, vars::{list, scalar})
   if [phi,vars]::listlist then
      Record(':-obj'= ':-CONSTRAINT', ':-expr'= phi, ':-vars'= vars)
   elif [phi,vars]::list(scalar) then
      thisproc([phi], [vars])
   else
      error "Invalid input. Arguments must be equipotent lists or both scalars."
   end if
end proc:

Also, a Record is a better data structure than a table for this purpose.

Since you want the value of x that gives the maximal value rather than the maximal value itself, you should use max[index], like this:

restart:
#Example functions:
(f,g):= 'unapply(randpoly([x,a], degree= 23, dense), [x,a])' $2;
a:= 1:
R:= Student:-Calculus1:-Roots(f(x,a), numeric);
R[max[index](g~(R,a))];

One must calculate coordinates, but that's easy enough. To avoid issues about the correspondence between sides and angles, I only handle the side-angle-side case in the procedure below. I assume that the angle is specified in degrees. Side1 is always the bottom, and Side2 is always the left side.

Triangle:= proc(
   Side1::positive, Side2::positive, 
   Angle::And(positive, satisfies(A-> A < 180))
)
   plot(
      [[0,0], [Side1,0], Side2*~[cos,sin](Angle*Pi/180), [0,0]],
      scaling= constrained, axes= none, _rest
   )
end proc:

As an example of using it
Triangle(2, 3, 65);

I immediately spot two anomalies in your code. The first is that the bound variable of maximize is s, and I don't see any s in the expression being maximized. This may or may not be the ultimate source of your error, but it needs to be corrected before we proceed.

The second anomaly isn't erroneous; it's just wasteful. You are forcing your Grid commands to be performed sequentially, not in parallel. Here's the relevant code:

i:= 1; while e < Grid[Map](errorterm, i) do i:= i+1 end do;

You are mapping over a single element. Each call to errorterm must finish before the next one begins. The above code is equivalent to

for i while e < errorterm(i) do end do;

In order to use Grid, you need to divide the algorithm into pieces that can be done in parallel, i.e., not sequentially. This is a fundamental theorethical limitation, not a limitation of Maple. So, if you're making multiple calls to Groebner where the input to one call doesn't depend on the output of another, then it can be done. But if you're making a single lengthy call to Groebner, then it can't.