Solution (source code)

= Solution

Iteration gives
$$
x_n=A^na+\sum_{j=0}^{n-1}A^jb
=A^na+(I-A^n)(I-A)^{-1}b,
$$
where the inverse exists by the hypothesis.