How to compute common zeros of system of polynomial equations with dimension 2?

edit

I have some ideal of homogenous polynomials defined over some finite field: J is the ideal of interest. However I can't call J.variety() since it is not zero dimensional. The system of equations might contain 6 or 231 polynomials in four variables:

  sage: [I.gens() for I in J.minimal_associated_primes()]
  [[x1 + x2, x0 + x3], [x2 + x3, x0 + x1], [x1 + x3, x0 + x2]]
  sage: J.dimension()
  2

Any ideas?

On request. The program computes the invariant under the group $2_{+}^{1+2\cdot2}$ homogenous polynomials of given degree. To construct the polynomials and get a list of them call: homInvar(6):

F = ZZ; 
a = matrix(F, [[0,0,0,-1],[0,0,1,0],[0,-1,0,0],[1,0,0,0]])
b = matrix(F, [[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]])
c = matrix(F, [[0,1,0,0],[-1,0,0,0],[0,0,0,1],[0,0,-1,0]])
d = matrix(F, [[0,1,0,0],[1,0,0,0],[0,0,0,-1],[0,0,-1,0]])

def getPart(deg):
    part = []
    for k in range(deg+1):
        for l in range(deg+1):
            for m in range(deg+1):
                for n in range(deg+1):
                    if k+l+m+n == deg:
                        part.append((k,l,n,m))
    return part

def p(x,part):
    x = x.list()
    pol = 0
    for p in [part]:
        mon = 1
        for i in range(len(p)):
            k = p[i]
            mon = mon * x[i]**k
        pol = pol + mon
    return pol

def getGroup():
    G = []
    for k in range(4):
        for l in range(4):
            for m in range(4):
                for n in range(4):
                    g = a**k*b**l*c**m*d**n
                    if G.count(g)==0:
                        G.append(g)
    return G

def reynolds(Gr,part):
    reyn = 0
    n = len(part)
    X = list(var('x%d' % i) for i in range(n))
    x = matrix([X]).transpose()
    for g in Gr:
        reyn +=  p(g*x,part=part)
        #print reyn
    return 1/len(Gr)* reyn

def homInvar(deg,Gr=getGroup(),getPart=getPart):
    parts = getPart(deg)
    inv = set([])
    for part in parts:
        r = reynolds(Gr,part)
        #print r
        if r != 0:
            inv.add(r)
    return inv