As used here, SymPy is just using Julia’s AbstractArray type. So this is a linear algebra question. Maybe, reduce(hcat, ([A^i for i in 0:(n-1)]..., B)) is what you want. (BTW, the code you had earlier seems to have an issue with your “D” matrix, so I couldn’t test.)