Here is an example for the finite field `GF(3)`

Let us first define the field:

```
sage: p = 3
sage: K = GF(p)
sage: K
Finite Field of size 3
```

Then, let us define the polynomial ring:

```
sage: R.<w,x,y,z> = K[]
sage: R
Multivariate Polynomial Ring in w, x, y, z over Finite Field of size 3
```

Then, let us define the 4 polynomials which have to vanish:

```
sage: eqn1 = 2*w*y^2 + y^2*z + z^3
sage: eqn2 = 4*x^3 - 2*x*y^2 - y^3
sage: eqn3 = 2*w^2*y - 2*x^2*y - 3*x*y^2 + 2*w*y*z
sage: eqn4 = w*y^2 + 3*w*z^2
```

Now, let us define the ideal generated by those polynomials:

```
sage: I = R.ideal([eqn1, eqn2, eqn3, eqn4])
sage: I
Ideal (-w*y^2 + y^2*z + z^3, x^3 + x*y^2 - y^3, -w^2*y + x^2*y - w*y*z, w*y^2) of Multivariate Polynomial Ring in w, x, y, z over Finite Field of size 3
```

Unfortunately, this ideal has dimension 1:

```
sage: I.dimension()
1
```

Hence, Sage is not able to find its variety (which corresponds to the set of common zeroes of the 4 polynomials above):

```
sage: I.variety()
ValueError: The dimension of the ideal is 1, but it should be 0
```

Hovewer there is a trick : on `GF(p)`

there fro the Frobenius, we know that the polynomial `x^p-x`

vanishes, hence, since we are looking to the solutions within `GF(p)^4`

, we can add those polynomials to the generators of the ideal:

```
sage: fw = w^p - w
sage: fx = x^p - x
sage: fy = y^p - y
sage: fz = z^p - z
sage: J = R.ideal([eqn1, eqn2, eqn3, eqn4, fw, fx, fx, fz])
sage: J.dimension()
0
sage: J.variety()
[{z: 0, y: 0, x: 0, w: 0}, {z: 0, y: 0, x: 0, w: 2}, {z: 0, y: 0, x: 0, w: 1}]
```

Do not use the symbolic ring SR, but a polynomial ring over a finite field.

To illustrate Frederic's answer:

BTW :

HTH,