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.

To do this with generality requires using diff rather than D, then stripping off the independent variable. Also, a procedure is required to declare names to be constants, so that they won't be interpreted as functions.

restart:
h:= F-> 
   subsindets(
      evalindets(
         diff(F(_),_)/F(_), 
         specfunc(diff), 
         d-> op(1,d)*'h'(op(1,d))
      ), 
      function(identical(_)), 
      f-> op(0,f)
    )
:
#Warning: DeclareConstant alters its argument, like `assume`.
DeclareConstant:= (C::evaln)-> assign(C, subs(__= C, _-> __)): 
DeclareConstant(a);
h((x+y)^a);                          

                             a (x h(x) + y h(y))
                             -------------------
                                    x + y       

Note that _ and __ are simply global names which should never be assigned a value.

This should help for the case of 2D plots. In your initialization file, include

plots:-setoptions(size= [1000,1000]);

adjusting the numbers to whatever is suitable to you.

Here's how you can get the solution using the Lagrange multiplier method in Maple:

restart:
V:= [X,Y]: #the variables
Vol:= Pi*X^2*(2*Y) + 2*Pi*X^2*(1-Y)/3:  #the objective function
Cons:= X^2+Y^2 - 1: #the constraint
#Solve for the Lagrangian being 0:
Sols:= [solve([map2(diff, Vol+lambda*Cons, V)[], Cons], {V[],lambda}, explicit)]:
Sols:= (S-> V=~S)~(map2(eval, V, Sols)):  #Strip off lambda.
#Select a critical point that maximizes the objective:
Sol:= Sols[max[index](eval~(Vol, Sols))]:
#Evaluate objective at that point:
<Sol[], Volume = expand(eval(Vol, Sol))>;

     X = (1/6)*sqrt(22+2*sqrt(13))
     Y = -1/6+(1/6)*sqrt(13)
     Volume = (35/81)*Pi+(13/81)*Pi*sqrt(13)

When is returns FAIL it means that Maple could not find a proof of the relation, nor could it find a counterexample. Your inequality is too complicated for Maple to affirm with absolute certainty. You'll need to try another approach. You may need to settle for a numeric verification, such as a plot.

Your command should be assume(2*y < x), not what you entered above.

There is a already a Maple command for counting the number of primes less than a given number: numtheory:-pi. It uses some more-advanced number theory than your procedure, so it's more efficient.

How is a table (using your terminology, not Maple's table) different from a Matrix? Note that a Maple Matrix can have any type of generic data as its entries.

MyTable:= Matrix((10,2), fill= 0);

I just included the fill= 0 to be analogous to the Mathematica; it's actually the default. If you want to get fancy, you can also include the embedded assignments to n and m:

`&=`:= proc(n::evaln, v) assign(n,v); v end proc:   
MyTable:= Matrix((n &= 10, m &= 2), fill= 0);

(I wonder why Maple's assign doesn't already return the value. That would be so useful, and would make it more compatible with other languages.)

[I began posting this before Markiyan's Reply appeared. Sorry for any duplication.]

The radius of the sphere is 6, so your plot is completely inside the sphere. To see the sphere, you need to make the r range larger. You also need to switch theta and phi: In Maple's spherical coordinates the second coordinate is the longitude, and the third coordinate is the latitude.

plots:-implicitplot3d(
   (r*sin(theta)*cos(phi))^2+(r*sin(theta)*sin(phi))^2+(r*cos(theta))^2=36,
   r=-7..7, theta=0..2*Pi, phi=0..Pi, coords= spherical
);

I know it's a hack, but you can plot the components separately. This is similar to Kitonum's, but uses a single plot command and avoids piecewise.

F:= x-> 
   if not x::realcons then 'procname'(args) 
   elif abs(x) < 2 then x
   else undefined
   end if
:
plot(
   [F~([seq](x-n, n= -8..8, 4))[], [seq]([[n-2,-2],[n+2,2]][], n= -8..8, 4)],  
   x= -10..10, scaling= constrained, color= blue, style= [line$5, point], 
   symbol= circle, symbolsize= 16
);

Since you've been here on MaplePrimes many times before and the solution to this system is straightforward (no singularities, convergence issues, etc.), I wonder what difficulty you encountered.

restart:
Eq1:= 
   diff(F(eta),eta$4) - 
   M*(eta*diff(F(eta),eta$3) + 3*diff(F(eta),eta$2) + 
      diff(F(eta),eta)*diff(F(eta),eta$2) - F(eta)*diff(F(eta),eta$3)
   ) - 
   Ha^2*diff(F(eta),eta$2)
:
Eq2:= 
   diff(G(eta),eta$2) + 
   Pr*M*(F(eta)*diff(G(eta),eta) - eta*diff(G(eta),eta)) +
   Pr*Ec*diff(F(eta),eta$2)^2 + Nb*diff(H(eta),eta) + 
   diff(G(eta),eta) + Nt*diff(G(eta),eta)^2
:
Eq3:= 
   diff(H(eta),eta$2) + M*Sc*(F(eta)*diff(H(eta),eta) - eta*diff(H(eta),eta)) +
   Nt*diff(F(eta), eta, eta)/Nb
: 
IC1:= F(0) = 0, (D@@2)(F)(0) = 0, D(G)(0) = 0, D(H)(0) = 0:
IC2:= F(1) = 1, D(F)(1) = 0, G(1) = 1, H(1) = 1:
 
params:= {Ec = .1, Nt = .1, Nb = .1, Sc = .5, Pr = 10, M = .5}:
HA:= [0, 2, 4, 6, 8]:

for k to nops(HA) do
   P||k:= plots:-odeplot( 
      dsolve(eval({Eq||(1..3), IC||(1..2)}, params union {Ha= HA[k]}), numeric),
      [[eta, F(eta)], [eta, G(eta)], [eta, H(eta)]],
      linestyle= [1,2,3],
      color= [red, blue, green, black, purple][k]
   )
end do:
plots:-display(P||(1..nops(HA)));

There's a simple trick that let's you exploit the internal representation of an Array to get what you want:

[rtable_elems(DataName)[]];

How about

plots:-polarplot([[r, 0, r= 0..1], [r, Pi/6, r= 0..1]], color= magenta, thickness= 6);

How about

plot([seq(-arctan(2*m*x)/(1-x^2), m= 1..10)], x= 0..10, discont);

It is always dangerous and unreliable to expect terms or factors to appear in a certain order. Isn't the result that you want the same as would would get from simply simplify(dairihensu1)?