Such 5-fold integrals containing discontinuous/nondifferentiable functions  are hard.
You must know to help Maple for them.
For example, use the method from my previous comment ==>

restart;
dd:= u=w*y..1, v=w*x..1, x=0..1, y=0..1, w=0..1:
int(abs((w*y-u)*(w*x-v)-w^2*x*y), dd);
evalf(%);

       -91/1080+4*ln(2)*(1/15)
       0.1005799889

Now, a numerical check is mandatory.

int(abs((w*y-u)*(w*x-v)-w^2*x*y)*piecewise(0 < -w*x*y+u*x and -w*x*y+u*x < 1, 1)*
piecewise(0 < -w*x*y+v*y and -w*x*y+v*y < 1, 1),
 [u = 0 .. 1, v = 0 .. 1, w = 0 .. 1, x = 0 .. 1, y = 0 .. 1],
 numeric, method=_CubaDivonne, epsilon=1e-5);

        0.100580642789240

We may conclude that the symbolic computation is correct [ 99.9% say :-) ].

restart:
with(Statistics):
#N :=3;# ... Some integer value >= 2;  
X := RandomVariable(Binomial(N, 1/2)): 
F:=CDF(X, s):
N:=3:
plot(F, s=-1..N+1, gridlines=true, axis[1]=[gridlines=N+1], discont, thickness=3);

A symmetric real matrix has orthogonal eigenvectors but the eigenvectors found by Eigenvectors are not necessarily orthogonal unless the eigenvalues are distinct or the matrix contains floats.
To obtain the orthogonal ones, use GramSchmidt.
In your case:

#  A = U . DiagonalMatrix(v) . U^+;

f := sin(x) / ( ( sin(a*x) )^a * ( sin((1-a)*x) )^(1-a) ):
MultiSeries:-limit(f, x=0,right) assuming a>0;

Same for a<0.

Maple is correct.

pde.mw

You have copied the expressions from the web page. Don't do this, because (unfortunately)  they are not designed to be copied.
Just open the help page (type ?Jacobian in Maple) and copy from there.
 

@Klausklabauter 

In this case all you need are the definitions and some Maple syntax.
Example

Maple sorts polynomials using monomial orders.
There is no such order with x^2 > y^2 > x*y.

Some tricks can be used e.g.

sort(x^2+y^2+x*y*``, order=plex(``,y),ascending);
  

but I don't like them because I prefer to use Maple for computations (and LaTeX for typesetting).

You should not expect 1.
Take the simpler A:=<1,2;3,-4>  to see why.

Maple computes the classical Cholesky decomposition A = L . L^*  (for a positive definite matrix A).  U is L^*.

It seems that you want the generalized one (the LDL decomposition)  A = LL . DD . LL^*

LL and DD can be obtained from L using:

S:=DiagonalMatrix(Diagonal(L));
DD:=S^2;
LL:=L.S^(-1);

Using limit for expressions containg floats is inherently problematic.
Compare:

limit(sin(1/3-x/3)/(3.*x-3.), x=1) ;
                         -0.1111111111
limit(sin(1/3.-x/3)/(3.*x-3.), x=1) ;
                         Float(undefined)

The  procedure  Basis  in  the  Groebner  package  computes  the  reduced Groebner  basis,  but after  this  computation  it  takes  the  primitive  parts  of all polynomials in  the basis. (Unfortunately this is not documented).

So, for
B := Basis(J,T) ;  
the "true"  reduced Groebner basis is
B/~LeadingCoefficient(B, T);

The integral cannot be computed symbolically.

restart;
F:=t -> Int(exp(-1-1/v)*(1-exp(v^2/(2*t^2)))/v^2, v = 0 .. 1):
evalf[30](F(1/10));
              -7.16224515512186261856578554808*10^18 
plot(F, 1..5);

evalf(Int(PDF(0.3*p1+0.7*p2, z, inert), z = 0 .. 1));
                          1.000000000
(only a few seconds).