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

To list the first 10 prime Leyland numbers, use

select(isprime, {seq(seq(a^b+b^a, a= 2..99), b= 2..99)})[1..10];

Some guess work was required to come up with the 99; it's much larger than is needed, but Maple's integer arithmetic and primality testing are so fast that it doesn't matter. Note that putting integers into a set automatically sorts them.

Here's a one-line solution:

cnt:= Statistics:-Tally(op~([g[]]), output= table);

This produces table output, like Joe's.