### Solving system of polynomial equation over finite field

Hello,

I am working in multivariate polynomial rings over small finite fields (with less then 15, or currently less than 10 elements). I have an algorithm, that outputs a system of equations. I need to find out if the system has a solution and if so, to find one. I am using calculation of groebner basis for this. But since my system typically has 20 or more variables and quite complicated polynomials, it seems that SINGULAR has problem to find the basis in reasonable time. Is there a way to speed up the calculation? For example, since I work in finite field, I can try all possible values for some variable. This simplifies the system and in some cases it seems to really help SINGULAR. But I can do that for only a very few variables.

Currently I work in PolynomialRing and before the grobner basis algorithm, I reduce my polynomials by FieldIdeal ~~(<x0**n -="" x0,="" ...,="" xm**n="" -="" xm="">) ~~(x_0^n - x_0, ..., x_m^n - x_m) of the polynomial ring and add this equations to the system, to force the solution in the finite field.
I tried to work in the quotient ring from the start, but it looks like it slowed down my whole algorithm and also slowed down the computation of the groebner basis (which is a little suprising for me).
I found that for GF(2) PolyBoRi is optimized for the kind of work I need, but I am not restricted to GF(2). Is there something similar for other finite fields?

I would appreciate any tips how to speed up the whole process. Thank you.