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.

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