Efficient solution to linear system involving Kronecker sum in MATLAB

High dimensional finite difference problems often lead to linear systems of the form $$A x = b, \qquad A = B_1 \oplus B_2 \oplus \cdots \oplus B_d,$$ where $$\oplus$$ denotes the Kronecker sum. $$B_i \in \mathbb{R}^{n \times n}$$ are 1D finite difference matrices, often tridiagonal or pentadiagonal. When $$d=2$$, for example, this is $$A = B_1 \otimes I + I \otimes B_2.$$

My question is about how to solve these systems most efficiently in MATLAB. The simplest approach is to construct the matrices using kron and solve with mldivide, but these systems can have a large bandwidth. In How to efficiently invert $K \otimes M+I_T\otimes \Sigma$? and in the literature on Sylvester equations, I've come across more efficient approaches based on diagonalizations or Schur decompositions of $$B_i$$. As $$B_i$$ are often Toeplitz or circulant, there are efficient FFT-based techniques for this, but I can't always make these structural assumptions. If my understanding is correct, the cost of this is $$O(d n^3 + n^d)$$.

Is there a more efficient way to solve these Kronecker sum systems when $$B_i$$ are banded? I speculate it is possible to do in $$O(d n + n^d)$$ time. Can this problem be decomposed into a sequence of 1D linear solves which are easily handled with mldivide?