For looking at the coefficients, you probably want to use

collect(_,'distributed');
coeffs(_,{y,z},_)

My notation might be off here. But collect distributed writes it as a sum of coefficients times monomials in y,z. Then coeffs will give you the coefficients, and store the corresponding monomials in an optional third argument.

Your 3 equations probably have no common solution for all paramaters. Any 2 of the equations should have a common solution (probably 4 of them). Probably solve will find them. You could solve on say the first two equations and the last two equations. This should give 4 solutions (each) in terms of your paramaters. You could then equate your parameters.

@_Maxim_ 

That's quite odd. In your example, if you remove the 'index' and replace one root with the negative you get the same issue as mine, i.e. duplicate (though correct, nonzero) eigenvectors.

Eigenvectors(<0, RootOf(_Z^2-2); -RootOf(_Z^2-2), 0>)

@Preben Alsholm 

I've found RootOf usually works fine without choosing which root. I believe the intended use is that any identical RootOf expression appearing in the same term is considered equal, such as how 'allvalues(m)' has 2 rather than 2^4 possible values.

The eigenvectors it outputs are correct in a way, they just don't correspond to the eigenvalues.

@Adam Ledger 

Congruence modulo polynomials is also handled in Maple. "rem" gives the congruence here - rem(a,b) will give a polynomial congruent to 'a' modulo 'b' of degree less than that of 'b'. "gcdex" computes the gcd of polynomials and can store the polynomials x and y such that xA+yB=gcd(A,B). If the gcd is 1, then x will be an inverse to A modulo B.

You could write functions poly_mult_mod(f,g,p) and poly_inverse_mod(f,p) using these. There should be a way of setting up infix notation if you desire.

There is a group package, which gives operations for the permutation groups for example.

@Kitonum 

Great thanks. vars is just a list of variables, it was intended to return, if vars=x, x[1]x[2]+x[3]x[4]+x[5]x[6]. A similar example with op() is

sum(op([[10],[10]][i]),i=1..2)

which returns 3. I see that

add(op([[10],[10]][i]),i=1..2)
sum('op'([[10],[10]][i]),i=1..2)

both return 20.

@tomleslie 

Initially, from copy pasting the output of showstat, the line "p1 := x^n+add(s || i*x^(n-i),i = 1 .. n);" gives a delimiters not matching error. Replacing it with cat(_,_) seems to work.

@Axel Vogt 
I want to know whether or not there are solutions over Q. Rank is not sufficient to determine solubility in my case, as there are more rows than columns (in A, for Ax=b), unless you mean something different?

@Carl Love 
Yes I'm looking for rational solutions. Solutions so far when they exist have had relatively small entries in numerator and denominator < 25. I have been using Matrix, does this effect the use of LinearSolve or just construction? The construction has been relatively quick, with LinearSolve being the bottleneck.

Do you know much about the modular method? Is it also quick at verifying when a solution does not exist?

Thanks for the tips, I'll try them out.

@Carl Love 
Cool thanks. Yes, the system is overdetermined. In any case I'm interested in, I've found that if a solution exists, then there is exactly one solution, I have only mild interest what this solution is, mostly I care whether or not a solution exists.

@Markiyan Hirnyk Thanks, but this is not helpful. The INVAR package mentioned only works for finite groups.

@Carl Love 

Thanks heaps, this is perfect.