All Maple's answers here are mathematically correct.
(It would be easier to see this for the similar function f(x) = Psi(1/x) - Psi(2/x)).

Facts:

• x = 0 is an essential non-isolated singularity (more exactly there is a sequence of poles at x_n, x_n < 0,  x_n --> 0) 
  So, the limit at 0 does not exist. The same (of course) for its derivative.
• The series given by MultiSeries is for x>0 (recall that it uses directional limits) and it is a "generalized series", not a Taylor one.
  (Note that the standard :-series refuses to give the expansion and invites to use asympt).

@Christopher2222 

No need to check, the expression is symmetric in x and y.

@Christopher2222 

A double integral must be computed:

with(Statistics):
X := RandomVariable(Normal(0, 1)): Y := RandomVariable(Uniform(-2, 2)):
#Probability(X*Y < 0);
f:=t->PDF(X,t);  g:=t->PDF(Y,t);
int(Heaviside(-s*t)*f(s)*g(t), s=-infinity..infinity, t=-infinity..infinity);
      1/2

@epostma 

But wouldn't this be just an ad hoc patch?
I think it would be interesting to investigate why those undefined and Dirac have appeared under int.
(And also why int is sometimes unable to compute an integral containing Heaviside without a convert/piecewise.)
 

Yes, I did not noticed that X3,X4 are uniform.
Then the CDF of Z is of course
CDF(X1, t)^2*CDF(X3, t)^2;

So, with assumptions Maple is correct.
Strange bug anyway.

@Christopher2222 

Unfortunately, all are wrong.
The correct CDF is:
CDF(X1,t)^4;
==>  (1/2+(1/2)*erf((1/2)*t*sqrt(2)))^4
(so, Mathematica's answer is ok).
Note that Maple also fails for  CDF(max(X1,X2,X3), t)
for which it results a limit at infinity < 0.6 (!)
but strangely, CDF(max(X1,X2), t) is correct.

@mmcdara 

- When seen as a distribution, the value of H at 0 does not matter.
- Before trying to compute HD, this entity must be defined as a distribution. How do you define it?

@mmcdara 

fsolve(CDF(X,x)=0.2, x);

@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.