DD[f_] := Expand[(f/.x->x+1)-f-(∂_y)f];
bas[n_] :=Table[(x^(n-k))(y^k), {k,0,n}];

bas[3]

DD/@bas[3]

(Subscript[MM, n_]):=Table[Coefficient[DD[a], b],{b,bas[n-1]},{a,bas[n]}]

(Subscript[MM, 3]) // MatrixForm

DD[1]