Use Gröbner bases to count the $3$-edge colorings of planar cubic graphs... I found a nice introduction on how to Use Gröbner bases to construct the colorings of a finite graph. 
Now my graphs $G=(V,E)$ are the line graphs planar of cubic graphs, so they are $3$-regular. The corresponding edge-adjacency matrices can be constructed, as shown here (in a crude way, I admit...). 
The existence of $3$-colorings on the edges of $G$ is quaranteed by planarity (On surfaces with higher genus only bipartite cubic graphs have chromatic index of $3$.)
Now let there be a field $F = \mathbb{Z}/3\mathbb{Z}$. Let's define two types of polynomials on $F$: 


*

*$f(x) = x(x-1)(x-2)= x^3-x=0$, which asks for one of three colors at edge $x_c$ and

*$g(x,y) = y^2+yx+x^2-1=0$, which asks for colors being different for the adjacent edges $y$ and $x$. 


Let $I$ be the ideal $I = (x_c^3-x_c \ |\ x_c \in E) + (x_r^2+x_rx_s+x_s^2-1 \ |\ x_{r,s} \in E)$. 
Note ${\displaystyle {\mathfrak {a}}+{\mathfrak {b}}}$ is the smallest left/right ideal containing both ${\displaystyle {\mathfrak {a}}}$ and ${\displaystyle {\mathfrak {b}}}$ (or the union ${\displaystyle {\mathfrak {a}}\cup {\mathfrak {b}}}$.

Moreover, every solution to this system yields a coloring and can be calculated by using the reduced Gröbner bases for the ideal $I \subseteq F[x_1, \ldots, x_{|E|}]$

Is it possible to calculate the number of solutions of the system of equations using Gröbner bases and if so how to do that?
 A: To count the number of solutions to a polynomial system over the algebraic closure of the field:
1) Compute a Groebner basis with respect to any ordering
2) For each variable $x_i$ there should exist a polynomial in the basis with leading monomial $x_i^{a_i}$.  If not, you have an infinite number of solutions.
3) You now need to count the number of monomials not divisible by a leading monomial from the Groebner basis.
For example, let $J := < x^2 + y + z, 2xy - z, xz - 5 >$;
We'll first compute a Groebner basis in lex order with $x > y > z$.  This is $\{z^4 - 10z^3 + 250, 10y - z^2, 50x + z^3 - 10z^2\}$.  The leading monomials are $\{z^4, y, x\}$ and the monomials which are not reducible are $\{1, z, z^2, z^3\}$.  The system has four solutions.  It works the same in any order.
In Maple:
with(PolynomialIdeals):
J := < x^2 + y + z, 2*x*y - z, x*z - 5 >;
NumberOfSolutions(J);

For your systems, you will want to add the option $(characteristic=3)$ when you construct the ideal.
Update #1
You could also use Singular from https://www.singular.uni-kl.de. Here is the example in Singular:
ring R = 0,(x,y,z),lp;
ideal J = [x^2 + y + z, 2*x*y - z, x*z - 5];
J = groebner(J);
kbase(J);

For your problem, you can change the ring definition to use characteristic 3 instead of 0 and replace "lp" (lex) with "dp" (grevlex).
