I found out that my computations boil down to the ‘’truncated’’ current Lie algebra $\mathfrak{gl}_N(\mathbb{C})\otimes\mathbb{C}[x]/(x^n)$, which is finite dimensional over $\mathbb{C}$, and the following code does what I need

```
d=3 #the size of matrices
n=12 #the ''depth'' of truncated polynomials, i.e. \mathbb{C}[x]/(x^n)
basis=''
for a in range(0,n):
for i in range(1,d+1):
for j in range(1,d+1):
basis=basis+',e'+'_'+str(i)+'_'+str(j)+'_'+str(a)
basis=basis[1:]
def rel(i,j,a,k,l,b):
if a+b>=n: return {}
elif k==j and i==l and i!=j: return {'e'+'_'+str(j)+'_'+str(j)+'_'+str(a+b):-1,'e'+'_'+str(i)+'_'+str(i)+'_'+str(a+b):1}
elif k==j and i==l and i==j: return {}
elif k==j and i!=l: return {'e'+'_'+str(i)+'_'+str(l)+'_'+str(a+b):1}
elif k!=j and i==l: return {'e'+'_'+str(k)+'_'+str(j)+'_'+str(a+b):-1}
elif k!=j and i!=l: return {}
relations={('e'+'_'+str(i)+'_'+str(j)+'_'+str(a),'e'+'_'+str(k)+'_'+str(l)+'_'+str(b)): rel(i,j,a,k,l,b) for i in range(1,d+1) for j in range(1,d+1) for k in range(1,d+1) for l in range(1,d+1) for a in range(0,n) for b in range(0,n)}
L = LieAlgebra(QQ, relations, names=basis)
ee=list(L.lie_algebra_generators())
def e(i,j,a):
if a<n: return ee[d*(i-1)+j+a*d^2-1]
else: return 0
```