from sympy import Matrix, Rational, mod_inverse, pprint
B = Matrix([[2, 2, 3, 2, 2], [2, 3, 1, 1, 4], [0, 0, 0, 1, 0], [4, 1, 2, 2, 3]]
    )
B_rref = B.rref(iszerofunc=lambda x: x % 5 == 0)
pprint(B_rref)
Sympy: Solving Matrices in a finite field