I'm trying to write a program that solves Sudokus using Gröbner bases. I introduced $81$ variables, $x_1$ to $x_{81}$, this is a linearisation of the Sudoku board.

The space of valid Sudokus is defined by:

for $i=1,\ldots,81$: $$F_i = (x_i - 1)(x_i - 2)\cdots(x_i - 9)$$ This represents the fact that all squares have integer values between 1 and 9.

for all $x_i$ and $x_j$ which are not equal but in the same row, column or block: $$G_{ij} = (F_i - F_j)/(x_i - x_j)$$ This represents that the variables $x_i$ and $x_j$ can not be equal.

All these $F_i$ and $G_{ij}$ together define the space of valid Sudokus. This consists of $891$ polynomials.

Now to solve a Sudoku we can add the clues to the space, so by example if the clue of a Sudoku is the first square is a $5$, then we add $(x_1 - 5)$ to the space. If we now take the Gröbner basis of this space we can directly see the solution for it.

I understand what I am doing this far. But I have trouble finding a computable manner for finding the gröbner bases. I have succesfully done everything for $4 \times 4$ Sudokus (or so-called Shidokus). But neither Maple nor Singular are giving me a result for the Gröbner basis of the $9 \times 9$ Sudoku space. You can see the commands I gave to Maple here. First I define the $891$ polynomials, then I ask for a basis of it.

I read papers saying it's feasible although non-performant to do what I strive for but I don't see how to find the solution, as they don't include many implementation details. Can anyone point me to a direction, making this problem easier for Maple or other software?