You can build the system matrix as follows A = RandomReal[{-1, 1}, {2, 2}]; A = A\[Transpose].A; n = 10; AA = Plus[ KroneckerProduct[SparseArray[{Band[{2, 1}] -> -1, Band[{1, 2}] -> 1}, {n, n}], A], KroneckerProduct[IdentityMatrix[n, SparseArray], N@IdentityMatrix[2]] ] The right-hand side can be assembled as follows: d = RandomReal[{-1, 1}, {n, 2}]; S1 = RandomReal[{-1, 1}, {2}]; Sn = RandomReal[{-1, 1}, {2}]; b = Flatten[d]; b[[;; 2]] += A.S1; b[[-2 ;;]] += A.Sn; Solve as usual (this is a usual linear system). I'd suggest a solver for banded matrices: LinearSolve[AA, b, Method -> "Banded"]