extend discretizer to work with matrix symbols
from discretizer import Discretizer
from discretizer import momentum_operators
from discretizer import coordinates
kx, ky, kz = momentum_operators
x, y, z = coordinates
A, B, C = sympy.symbols('A B C', commutative=False)
sigma = sympy.MatrixSymbol('Sigma', 2, 2)
hamiltonian = sympy.Matrix([[kx * A(x) * kx, B(x,y)*kx], [kx*B(x,y), C*ky**2]]) + sigma
tb = Discretizer(hamiltonian, discrete_coordinates={'x', 'y'}, lattice_constant=2.0, verbose=True)