Here is what comes first in mind:
sage: R.<a,b,c,d,e,f> = GF(2)[]
sage: I = R.ideal([1+a+c+e, b+d+f, 1+a*e+c*e+a*c, b*e+a*f+c*f+e*d+b*c+a*d, b*f+f*d+b*d])
sage: I.variety()
ValueError: The dimension of the ideal is 2, but it should be 0
Unfortunately, even if we are in finite fields, Sage can only enumerate solutions of a system of dimension zero.
Let us, for example add the two equations a=1 and b=1:
sage: I = R.ideal([1+a+c+e, b+d+f, 1+a*e+c*e+a*c, b*e+a*f+c*f+e*d+b*c+a*d, b*f+f*d+b*d, a+1, b+1])
sage: I.variety()
[]
Which means that there are no solution in GF(2), but, as you mentionned, there are solutions in GF(4):
sage: I.variety(GF(4,'w'))
[{b: 1, c: 1, a: 1, f: w, d: w + 1, e: 1},
{b: 1, c: 1, a: 1, f: w + 1, d: w, e: 1}]
If you want all solutions, instead of arbitrary equations a=1, you could add equations which are tautological in GF(4), namely, the equation given by the Frobenius:
sage: [x^4+x for x in R.gens()]
[a^4 + a, b^4 + b, c^4 + c, d^4 + d, e^4 + e, f^4 + f]
So, you can define:
sage: I = R.ideal([1+a+c+e, b+d+f, 1+a*e+c*e+a*c, b*e+a*f+c*f+e*d+b*c+a*d, b*f+f*d+b*d] + [x^4+x for x in R.gens()])
sage: I.variety(GF(4,'w'))
[{b: 0, c: w, a: w + 1, f: 0, d: 0, e: 0},
{b: 0, c: w + 1, a: w, f: 0, d: 0, e: 0},
{b: 0, c: 1, a: 1, f: 0, d: 0, e: 1},
{b: 0, c: 0, a: w + 1, f: 0, d: 0, e: w},
{b: 0, c: w + 1, a: 0, f: 0, d: 0, e: w},
{b: 0, c: 0, a: w, f: 0, d: 0, e: w + 1},
{b: 0, c: w, a: 0, f: 0, d: 0, e: w + 1},
{b: w + 1, c: w + 1, a: w, f: 1, d: w, e: 0},
{b: w, c: w, a: w + 1, f: 1, d: w + 1, e: 0},
{b: w + 1, c: 1, a: 1, f: 1, d: w, e: 1},
{b: w, c: 1, a: 1, f: 1, d: w + 1, e: 1},
{b: w + 1, c: 0, a: w + 1, f: 1, d: w, e: w},
{b: w, c: w + 1, a: 0, f: 1, d: w + 1, e: w},
{b: w + 1, c: w, a: 0, f: 1, d: w, e: w + 1},
{b: w, c: 0, a: w, f: 1, d: w + 1, e: w + 1},
{b: w + 1, c: w, a: w + 1, f: w, d: 1, e: 0},
{b: 1, c: w + 1, a: w, f: w, d: w + 1, e: 0},
{b: w + 1, c: 1, a: 1, f: w, d: 1, e: 1},
{b: 1, c: 1, a: 1, f: w, d: w + 1, e: 1},
{b: w + 1, c: w + 1, a: 0, f: w, d: 1, e: w},
{b: 1, c: 0, a: w + 1, f: w, d: w + 1, e: w},
{b: w + 1, c: 0, a: w, f: w, d: 1, e: w + 1},
{b: 1, c: w, a: 0, f: w, d: w + 1, e: w + 1},
{b: w, c: w + 1, a: w, f: w + 1, d: 1, e: 0},
{b: 1, c: w, a: w + 1, f: w + 1, d: w, e: 0},
{b: w, c: 1, a: 1, f: w + 1, d: 1, e: 1},
{b: 1, c: 1, a: 1, f: w + 1, d: w, e: 1},
{b: w, c: 0, a: w + 1, f: w + 1, d: 1, e: w},
{b: 1, c: w + 1, a: 0, f: w + 1, d: w, e: w},
{b: w, c: w, a: 0, f: w + 1, d: 1, e: w + 1},
{b: 1, c: 0, a: w, f: w + 1, d: w, e: w + 1}]
This can be done in any finite extension of GF(2) using the same trick, for example:
sage: n = 2^6
sage: I = R.ideal([1+a+c+e, b+d+f, 1+a*e+c*e+a*c, b*e+a*f+c*f+e*d+b*c+a*d, b*f+f*d+b*d] + [x^n+x for x in R.gens()])
sage: I.variety(GF(n,'w'))
[here you get a dictionary of 8191 solutions]
Now, if you want to go to the algebraic closure, you can not solve the dimension issue as before, so, if you have enough equations (e.g adding a+c and b^8+b), you can do:
sage: K = GF(2).algebraic_closure()
sage: I = R.ideal([1+a+c+e, b+d+f, 1+a*e+c*e+a*c, b*e+a*f+c*f+e*d+b*c+a*d, b*f+f*d+b*d, a+c, b^8+b])
sage: I.variety(K)
[{b: 0, c: 1, a: 1, f: 0, d: 0, e: 1},
{b: z3^2 + z3, c: 1, a: 1, f: z6^2 + z6, d: z6^5 + z6^2, e: 1},
{b: z3^2 + z3 + 1, c: 1, a: 1, f: z6^3, d: z6^5 + z6^3 + z6 + 1, e: 1},
{b: 1, c: 1, a: 1, f: z2, d: z2 + 1, e: 1},
{b: 1, c: 1, a: 1, f: z2 + 1, d: z2, e: 1},
{b: z3, c: 1, a: 1, f: z6^4 + z6^2, d: z6^5 + 1, e: 1},
{b: z3 + 1, c: 1, a: 1, f: z6^4 + z6^3 + z6 + 1, d: z6^5 + z6^3 + z6^2 + z6 + 1, e: 1},
{b: z3, c: 1, a: 1, f: z6^5 + 1, d: z6^4 + z6^2, e: 1},
{b: z3^2 + z3, c: 1, a: 1, f: z6^5 + z6^2, d: z6^2 + z6, e: 1},
{b: z3^2, c: 1, a: 1, f: z6^5 + z6^2 + z6 + 1, d: z6^5 + z6^4, e: 1},
{b: z3^2 + 1, c: 1, a: 1, f: z6^5 + z6^3 + 1, d: z6^5 + z6^4 + z6^3 + z6^2 + z6 + 1, e: 1},
{b: z3^2 + z3 + 1, c: 1, a: 1, f: z6^5 + z6^3 + z6 + 1, d: z6^3, e: 1},
{b: z3 + 1, c: 1, a: 1, f: z6^5 + z6^3 + z6^2 + z6 + 1, d: z6^4 + z6^3 + z6 + 1, e: 1},
{b: z3^2, c: 1, a: 1, f: z6^5 + z6^4, d: z6^5 + z6^2 + z6 + 1, e: 1},
{b: z3^2 + 1, c: 1, a: 1, f: z6^5 + z6^4 + z6^3 + z6^2 + z6 + 1, d: z6^5 + z6^3 + 1, e: 1}]
