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.

Did you derive the polylog expression from an alternating series? If so, then you can use the fact that the absolute value of the approximation error is less than the absolute value of the first term of the series that you don't use. Just add up the first however many terms of the series that you need.

Threads shares the memory among the several processes. Thus global variables, if used at all, must be used very carefully, lest multiple processes update them in the wrong order. Most of the Maple system was written long before the Threads package, and the use of globals (often in the form of "environment variables") is too deeply engrained in the system for it to be changed now. The help page ?index,threadsafe lists the Maple commands that can be used with Threads. That list does not include int, sum, plot, floor, or piecewise---all of which you use. When you use commands that aren't "thread safe", the results are exactly as you described: random, unpredictable, and with frequent kernel crashes.

There is another multiprocessing package that doesn't use shared memory: Grid. The overhead for using it is much higher than with Threads, so the payoff is less, and of course the memory usage is higher.

The volume is the integral of Pi*f(x)^2 over the x interval. For example:

f:= x-> sqrt(1-x^2); #semicircle
Int(Pi*f(x)^2, x= -1..1);
value(%);

The arclength of f(x) on the interval x= a..b is the integral of sqrt(1+diff(f(x),x)^2) over a..b. For example:

f:= x-> sqrt(1-x^2); #semicircle
Int(sqrt(1+D(f)(x)^2), x= -1..1);
value(%);

Download SurfaceInt.mw

The command is RootFinding:-Analytic.

RootFinding:-Analytic(
   u*(BesselJ(0,u)^2 + BesselJ(1,u)^2) - 2*BesselJ(0,u)*BesselJ(1,u),
   u, re= -99..99, im= -99..99
);

The limitation is that you need to specify a rectangle over which to find the roots.

Using Maple 2016. I have no problem with your code. I don't have Maple 2015 to test. However, the with(LinearAlgebra) is unnecessary, so try removing it, just in case.

Regarding your question about solve applied to a vector: It takes a lot less than a few extra lines of code to do this. All that's needed is

solve(convert(v, set), {x,y});

How about

FallingFactorial:= (x,n)-> pochhammer(x-n+1, n);