@JAMET 

u runs in the set {0, 1/n, ..., a/b, ..., 1}  constructed just above it.

Please explain for this simpler expression.

f:=diff(v(x, y), x, y) * diff(u(x, y), y, y)^2;

What are you expecting for int(f, x=0..1, y=0..1); ?

@Carl Love 

Are you saying that the global variables are saved in a .m file without being explicitely present in the arguments of save?

I have not seen this behavior. Probably you mean something else.

@acer 

Thank you, good to know. Probably I should have guessed it. Actually, if B is not definite, some eigenvalues could be undefined (even if symmetric and commuting A,B implies real eigenvalues -- but their number could be < n).

@Carl Love 

You are right, I have assumed (without checking) that A,B commute (A.B=B.A).

Edit. Even when A,B commute and have hfloats, the eigenvalues still have imaginary 0s.
So, the algorithm in this case ignores the symmetry.

@Carl Love 

Yes, M is symmetric, but it seems that this is ignored by Eigenvalues, because 0 imaginary parts still appear. Probably the symmetry is taken into account only for hfloats.

@tomleslie

Download JK.mw

@tomleslie 

The OP does not want to reassign anything; he wants some computations with U.

@mmcdara 

You missed 13 cycles.

{[1,3,5],[1,3,6],[1,4,5],[1,4,6],[2,3,5],[2,3,6],[2,4,5],[2,4,6],[1,3,2,4],[1,3,2,5],[1,3,2,6],[1,3,5,4],[1,3,6,4],[1,4,2,5],[1,4,2,6],[1,5,2,6],[1,5,3,6],[1,5,4,6],[2,3,5,4],[2,3,6,4],[2,5,3,6],[2,5,4,6],[3,5,4,6],[1,3,2,4,5],[1,3,2,4,6],[1,3,2,5,4],[1,3,2,6,4],[1,3,5,2,4],[1,3,5,2,6],[1,3,5,4,6],[1,3,6,2,4],[1,3,6,2,5],[1,3,6,4,5],[1,4,2,3,5],[1,4,2,3,6],[1,4,5,2,6],[1,4,5,3,6],[1,4,6,2,5],[1,4,6,3,5],[1,5,2,3,6],[1,5,2,4,6],[1,5,3,2,6],[1,5,4,2,6],[2,3,5,4,6],[2,3,6,4,5],[2,4,5,3,6],[2,4,6,3,5],[1,3,2,5,4,6],[1,3,2,6,4,5],[1,3,5,2,4,6],[1,3,5,2,6,4],[1,3,5,4,2,6],[1,3,6,2,4,5],[1,3,6,2,5,4],[1,3,6,4,2,5],[1,4,2,5,3,6],[1,4,2,6,3,5],[1,4,5,2,3,6],[1,4,5,3,2,6],[1,4,6,2,3,5],[1,4,6,3,2,5],[1,5,3,2,4,6],[1,5,4,2,3,6]}

@weidade37211 

When working with such "extreme" numbers a good idea is to transform the equation symbolically.
 

@Carl Love 

Here is the winding number version. It should work for any closed curve.

ContourInt:= proc(
   f::algebraic, 
   C::(name= algebraic), 
   trange::name= range(realcons)
)
local z:= lhs(C), t:= lhs(trange), r:= rhs(C), eq, R;
R:=remove(z-> type(rhs(z), infinity), op~({singular(f,z)} ));
2*Pi*I*add(round(evalf(Int(diff(r,t)/(r-eval(z,eq)),trange)/(2*Pi*I)))*residue(f,eq),eq=R);
end proc:

@Carl Love 

To decide whether a singular point is in the interior we can compute its winding number (it's enough to approximate and apply round).

@Lali_miani 

You already got the explanation: sqrt(2*x+1)/(2*x+1) is automatically simplified to 1/sqrt(2*x+1). The automatic simplifications cannot be prevented, e.g. you cannot see 2/10 instead of 1/5. There are some tricks, but I think that it is much more important to focus on maths!

@Carl Love 

This is probably the most efficient method.
But I think that for a beginner, a more useful construction (to be often used) is

MyMatrix := Matrix( 2,3, (i,j) -> elem[i,j](t) );

@dharr 

So, I'd formulate:

When x is a sparse Matrix, Vector or rtable, for efficiency reasons, only the nonzero entries are scanned and the zero entries are then appended.