How to solve system of equations in polynomial ring over GF with a different variables

edit

solve() is a symbolic solve function, which works over the symbolic ring SR not polynomial rings. Since you work with polynomials, it's better to employ the polynomial ideal machinery instead:

F = GF(31)
R.<x, y, z> = PolynomialRing(F)
system = [
    x - (3 * y),
    z - (2 * x)
]
J = Ideal(system)
print( J.variety() )

However, in this case you'd get an error:

ValueError: The dimension of the ideal is 1, but it should be 0

which tells us about the system having too many "degrees of freedom". Since we work in a finite field, it is easy to deal with it by fixing one of the variables, say, x:

sols = sum( (Ideal(system+[x-x0]).variety() for x0 in F), [] )
print(sols)

which prints the solutions:

[{z: 0, y: 0, x: 0}, {z: 2, y: 21, x: 1}, {z: 4, y: 11, x: 2}, {z: 6, y: 1, x: 3}, {z: 8, y: 22, x: 4}, {z: 10, y: 12, x: 5}, {z: 12, y: 2, x: 6}, {z: 14, y: 23, x: 7}, {z: 16, y: 13, x: 8}, {z: 18, y: 3, x: 9}, {z: 20, y: 24, x: 10}, {z: 22, y: 14, x: 11}, {z: 24, y: 4, x: 12}, {z: 26, y: 25, x: 13}, {z: 28, y: 15, x: 14}, {z: 30, y: 5, x: 15}, {z: 1, y: 26, x: 16}, {z: 3, y: 16, x: 17}, {z: 5, y: 6, x: 18}, {z: 7, y: 27, x: 19}, {z: 9, y: 17, x: 20}, {z: 11, y: 7, x: 21}, {z: 13, y: 28, x: 22}, {z: 15, y: 18, x: 23}, {z: 17, y: 8, x: 24}, {z: 19, y: 29, x: 25}, {z: 21, y: 19, x: 26}, {z: 23, y: 9, x: 27}, {z: 25, y: 30, x: 28}, {z: 27, y: 20, x: 29}, {z: 29, y: 10, x: 30}]

ADDED. As John pointed out, solutions of the form like x = 3 * y, z = 6 * y can be seen from the generators of the ideal:

print(Ideal(system).gens())

which gives

[x - 3*y, -2*x + z]

that is, y = x/2 and z = 2*x.

The solve_mod command might be useful:

Return all solutions to an equation or list of equations modulo the given integer modulus. Each equation must involve only polynomials in 1 or many variables.

Type solve_mod? for the help message and examples. This is most useful when there are relatively few solutions, maybe even when there are relatively few possible solutions. It is very slow for the particular equations you provide, partly because the modulus is 31 and it may very well be trying all possible values for x, y, and z. Also, it will list all solutions, not solve it symbolically the way you are asking.