@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

It's a pity that fsolve fails for the equation erf(x) = 1 - 10^(-20)  when Digits = 5000.
It works when Digits = 4000.
The equation can be solved via

Digits:=5000:
x1:=RootFinding:-Analytic(erf(x) = 1 - 10^(-20), 6-I/10 .. 7+I/10);

or even directly with the good old Newton's method:

restart;
y1:=1-10^(-20):
Digits:=5100:
a:=6.;h:=1. : N:=0:
K:=evalf(sqrt(Pi)/2):
while  abs(h)>1e-5000 do
h:=(erf(a)-y1)*exp(a^2)*K;
a:=a-h:  N:=N+1;
od:  N; a;