@Markiyan Hirnyk 

My comment was towards the rest of the Maple community because I have anticipated the (usual) title of your reply.
The past demonstrates that our maths, Maple and communication manner are too different. 

The Multiseries:-limit seems to differ from the usual limit in some aspects.

"The variable is assumed to tend to its limit along a ray coming from the origin."

The "directional" limits can be strange. There exist analytic functions in a punctured disc where all directional limits exist but the (global) limit does not.

MultiSeries:-limit(sin(x)/x, x = infinity, complex);
       0

Mathematically should be undefined. But if infinity is the real +oo and the limit is directional then it could be interpreted as correct (but then the complex option is simply superfluous). [Interpreted this way, both OP's limits are correct].

Note that

MultiSeries:-limit(sin(x)/x, x = infinity+I*infinity, complex);
returns unevaluated.

In this case the usual limit returns unevaluated
which means "I don't know"  for x = infinity
but gives an error (Error, (in limit) invalid limiting point)
for x = infinity + I*infinity.

Edit.
My assertion about "directional" (radial) limit in MultiSeries and the complex option is confirmed by the fact that
MultiSeries:-limit((x^2-1)*sin(1/(x-1)), x = -infinity, complex);
returns  - infinity.
MultiSeries:-limit((x^2-1)*sin(1/(x-1)), x = infinity*I, complex);
returns   infinity*I.

while

limit((x^2-1)*sin(1/(x-1)), x = -infinity, complex);
returns  -infinity+infinity*I

@acer 

The problem seems to have historical roots. There was a time when Maple did not have strings, but only names (symbols).
For keywords it is of course easier to write (and read) coords = polar than "coords" = "polar".
But if we are going to be forced to use
':-coords' = ':-polar'
(because the number of used modules increases)
then a good idea would be to allow everywhere strings.
 

@Carl Love 

On my computer the speed increases only by 5%.
But the real problem is that you computed a wrong matrix
(because n <> 2*n+1).

Here is a version where the speed increases 10 times!

n:= 1000:
A:= Matrix((2*n+1 $ 2), datatype= float[8]):
V:=Vector(2*n+1,datatype= float[8]):
evalhf(
   proc(A,V,n)
   local 
      i, j, u, h:= 5/n;
      for i to 2*n+1 do u:=(i-1-n)*h; V[i]:=`if`(u=0,1.0, (sin(u)/u)^2) od;
      for i to 2*n+1 do
         for j from i to 2*n+1 do
            A[i,j]:= V[i]*V[j]
         end do
      end do
   end proc
   (A,V,n)
):
A:= Matrix(A, shape= symmetric):

Probably a compilation would further speed up the computations.

 

@tomleslie 
With a couple of little changes the timing is almost the same:

sinc2:= proc(x) option remember; (sin(x)/x)^2 end:
sinc2(0):=1: sinc2(0.):=1.:
n:=1000;
a:=5.0; h:=a/n;
Z:=Matrix(2*n+1, (i,j) -> sinc2((i-1-n)*h)*sinc2((j-1-n)*h), datatype=float[8],shape=symmetric);

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

The next one is 30% faster:

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

@Markiyan Hirnyk 

Cite:

Description
The implicitplot command computes the two-dimensional plot of an implicitly defined curve.  By default, the curve is computed in Cartesian coordinates.

Note that I have added "probably" and "experimental".
That's the "reality".
 

In the help it is mentioned that implicitplot is designed for curves.
It's true that in an example appears a region (f(x,y)<0) but probably this part is only experimental.

@jessica99 

truncation order.

You could have typed ?mtaylor  for Maple help.

Anyway implicitplot should be used only when explicit or parametric are out of the question. So:

plot([seq(arctan(surd(3*cos(x),3))+k*Pi,k=-1..1)], x=-Pi..Pi,y = -(1/2)*Pi-1 .. (1/2)*Pi+1, color=red) ;

@hitstudent 

Maple has it but the user must formulate mathematically the problem.
In this problem Maple will by really helpful if e.g. the direction of the projection is not the normal of the plane or if the sides of the square are not paralel/orthogonal to the base of the cylinder.

@Markiyan Hirnyk 

Thank you. Not a bad job done by Mathematica.
However the two samples N(1,1),  N(3,1/2) seem to be too separated, so the task is easier.

@Carl Love 

I don't think that the command works this way. It could return a convex combination of a normal and some other distribution. It would be interesting to merge two equal samples N(m1,s1), N(m2,s2) and see whether m1,m2,s1,s2 are retrived.

Maybe Markiyan could do this.

@Markiyan Hirnyk 

It was just an opinion, I am not a statistician. However, it does not seem to be  a very complicated command. In your example it returned a convex combination of 2 standard distributions. Having a finite set F of such distributions, it is not difficult to look for a distribution c*f + (1-c)*g,  f,g in F, c in [0,1]  and determine the pair (including parameters) which best fits for a given sample.

@TheFixer 

For cones and cylinder "h" is not the same. Look at my vol.

@TheFixer 

In my formula the max volume is 3.175419717

while X=Y=sqrt(2)/2 corresponds to t=Pi/4 and the volume is 2.528158531.