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);

Change the command to

plots:-fieldplot(<x,y>, x=0..1, y=0..1, fieldstrength= ':-log');

This'll work in either case, meaning regardless of whether RealDomain has been loaded.

With any keyword option to any command, the keyword keyword can always be replaced by ':-keyword', and the code will be more robust if you do so. The :- means to use the global version of the name. In this case, that means regular log instead of RealDomain:-log. The ' ' means to ignore any value that may have been assigned to the name. It's not necessary in this particular situation, but it's a good idea in general when you're using these keywords inside a procedure that may not be aware of its global environment.

This is solved by Maple just like any other symbolic BVP:

Sol:= dsolve({ode, bcs}):
simplify(Sol) assuming x > y;
simplify(Sol) assuming x <= y;

I find Kitonum's Answer using remove a bit too ad hoc. And I'd prefer that you not use `error` (necessarily with the quotes!), although technically it works, so I'll use Error in my example:

simplify(Error, {h^6 = 0});

This should annihilate h^6 terms and also any terms with higher powers of h.

VV wrote:

• Note that evalhf does not change much because Z has float[8].

With the correct use of evalhf, your best time can be reduced by another 27%.

restart:
st:= time():
n:= 1000:
a:= 5.0: h:= a/n:
sinc3:= proc(i) option remember; local u:=(i-1-n)*h; (sin(u)/u)^2 end:
sinc3(n+1):= 1:
Z:=Matrix(2*n+1, (i,j) -> sinc3(i)*sinc3(j), datatype=float[8],shape=symmetric):
time()-st;
                             6.359

restart:
st:= time():
n:= 1000:
A:= Matrix((2*n+1 $ 2), datatype= float[8]):
evalhf(
   proc(A,n)
   local 
      i, j, u, h:= 5/n, N:= 2*n+1,
      sinc3:= (i,n,h)-> (u-> `if`(u=0, 1, (sin(u)/u)^2))((i-1-n)*h)
   ;
      for i to N do
         for j from i to N do
            A[i,j]:= sinc3(i,n,h)*sinc3(j,n,h)
         end do
      end do
   end proc
   (A,n)
):
A:= Matrix(A, shape= symmetric):
time()-st;         
                             4.625

(6.390-4.625)/6.390;
                       0.276212832550861

My guess is that you want the Maclaurin series as a symbolic summation. In that case use

convert(erf(x), FormalPowerSeries);

The clock resolution is as you say. Precisely, it is 1/64 s. I've occasionally seen systems where it's 1/256 s, but never finer than that. To compensate for this, you need to call the code repeatedly, preferably with randomized input, then divide by the number of trials. There are several other subtleties involved in accurate timing, such as accounting for the garbage collection, the memory allocation, and the initial reading of code from libraries. I've written a great many articles on MaplePrimes (and earlier Maple-related fora) on this over the years. I'll give more details later.

Someone may suggest using CodeTools:-Usage with the iterations option. This may be adequate for some very simple pieces of kernel code. The lack of randomized input should make any experimental scientist skeptical about the generalizability of any conclusions.

To answer your question #2 only: The methods are listed at ?dsolve,numeric. They are rkf45, ck45, rosenbrock, bvp, rkf45_dae, ck45_dae, rosenbrock_dae, dverk78, lsode, gear, taylorseries, mebdfi, and classical. These all have individual help pages, which are hyperlinked from ?dsolve,numeric. Many of these have submethods; for example, euler is a submethod of classical.