//Demo modeling example

q := 2;
n := 5;
K := FiniteField(q);

R := BooleanPolynomialRing(2*n^2, "grevlex");

vars := ["a_" * IntegerToString(i) * "_" * IntegerToString(j) : j in [1..n], i in [1..n] ];
vars := vars cat ["b_" * IntegerToString(i) * "_" * IntegerToString(j) : j in [1..n], i in [1..n] ];
AssignNames(~R, vars);

A_ := Matrix(R, n, n, [R.((i-1)*n+j) : j in [1..n], i in [1..n]]);
B_ := Matrix(R, n, n, [R.((i-1)*n+j+n^2) : j in [1..n], i in [1..n]]);



C1 := RandomMatrix(K,n,n);
while Rank(C1) lt 4 do
    C1 := RandomMatrix(K,n,n);
end while;
C2 := RandomMatrix(K,n,n);
while Rank(C2) lt 4 do
    C2 := RandomMatrix(K,n,n);
end while;
A := Random(GL(n, K));
B := Random(GL(n, K));
D1 := A*C1*B;
D2 := A*C2*B;

D1_var := A_*C1*B_;
D2_var := A_*C2*B_;

system := [];
for i in [1..n] do
    for j in [1..n] do
        Append(~system, D1[i][j]-D1_var[i][j]);
        Append(~system, D2[i][j]-D2_var[i][j]);
    end for;
end for;

//try to solve system
SetVerbose("Groebner", 0);
I := Ideal(system);
GroebnerBasis(I);
Variety(I);