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

A bug in sum  occurs for the inner sum.

x:=3:
sum(binomial(x,j)*(1/3)^j*(2/3)^(x-j), j = i+1 .. x);
       1

Even for a simpler sum:

sum(binomial(10,j), j = k .. 10);
                            1024

Using infolevel[sum]:=5;  it appears that the bug is in the summation using hypergeometric functions.

Just replace rotate with 'rotate'.

Your n-dimensional integral

J:= n -> int(add(x[i],i=1..n)^2, [seq(x[i]=0..1,i=1..n)]);

has the exact value

Jexact:= n -> n*(3*n + 1)/12;

seq('J'(n)=Jexact(n), n=1..20);
J(1) = 1/3, J(2) = 7/6, J(3) = 5/2, J(4) = 13/3, J(5) = 20/3, J(6) = 19/2, J(7) = 77/6, J(8) = 50/3, J(9) = 21, J(10) = 155/6, J(11) = 187/6, J(12) = 37, J(13) = 130/3, J(14) = 301/6, J(15) = 115/2, J(16) = 196/3, J(17) = 221/3, J(18) = 165/2, J(19) = 551/6, J(20) = 305/3

(Just in case you want to test the Monte Carlo method for a higher n.)

restart;
J:=Int((-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2), x = 4*(1/10) .. 6*(1/10)):
convert(J,Sum):
evalf[15](%);

                      -0.00555168769360751