Yes, it is somewhat trickier. For a matrix of floats, roundoff error usually turns an eigenvalue with multiplicity > 1 into several close, but not exactly equal, eigenvalues. Exactly how this will happen, and what will be the resulting eigenvectors, is probably not easily predictable.

Any nonzero constant multiple of an eigenvector is also an eigenvector. The real question, I think, is "how does Maple normalize its eigenvectors?" From a few tests (using LinearAlgebra[Eigenvectors]), it seems to me (for a Matrix with floats, possibly complex) that Maple chooses an eigenvector of length 1 (i.e. the sum of the squares of the absolute values of the entries is 1), and the entry that is largest in absolute value (or one of those in case of a tie) is real. That entry can be either positive or negative, so this isn't a complete specification, but perhaps it's a start.

Perhaps that's hinted to by the following cryptic statement in the ?rtable_elems help page: Only the elements that explicitly need to be specified, taking into account storage, are returned. In any case, it does seem that (even though storage is rectangular rather than sparse for this Array) zeros are not returned by rtable_elems, and in fact the examples in the help page show this. I don't know if there's a work-around [other than kludgy things like adding epsilon to all the elements first].

... or in John's application,

> map([lhs], select(t -> type(rhs(t),posint), rtable_elems(A)));

Yes, that's a better way of doing it. Somehow I missed rtable_elems.

Your "christoffel" is not an Array. It is a module. AFAIK there is no problem in converting a 3-dimensional Array to an array.

How about this way?

> select(t -> type(A[op(t)], posint),   {indices(convert(A,array))});

It works fine for me in Maple 9.5. Can you show us a specific example of this problem?

It works fine for me in Maple 9.5. Can you show us a specific example of this problem?

Or even simpler: [seq(seq(seq([a,b,c],a=A),b=B),c=C)];

Or even simpler: [seq(seq(seq([a,b,c],a=A),b=B),c=C)];

That problem with < bites again... I think what you meant was something like

> inequalities({x^2 + y^2 <= 2, x^2 + y^2 >= 1},
    x = -2 .. 2, y = -2 .. 2);

That problem with < bites again... I think what you meant was something like

> inequalities({x^2 + y^2 <= 2, x^2 + y^2 >= 1},
    x = -2 .. 2, y = -2 .. 2);

Allan, I think what you're seeing there for small values of r is just roundoff error, exacerbated by the fact that floats are being used. Consider: > rxpr:= convert(xpr, rational); series(rxpr, r); # returns quickly > evalf(%); series(2.720294101+.2646772639e-1*r^2-.5606773957*r^4+O(r^5),r,5) > limit(rxpr, r=0); # blows up

Perhaps you're looking for something like this? > Data := [[1,3],[2,1],[3,2],[4,3],[5,5]]; P := t -> a + b*t + c*t^2; L := map(u -> u[2] - P(u[1]), Data); Res := Optimization:-LSSolve(L); plots[pointplot](zip((u,r) -> [u[1],r], Data, subs(Res[2], L)));

Of course it's not the same thing. The rules for 2D math entry in Standard GUI are different. That's how you can get all those fancy characters. This is the way it was designed to work.