Author here. It actually is a block tridiagonal matrix with tridiagonal blocks. The test problem only has a quarter millions of unknowns so, for a production run, I wouldn't bother write my own sparse solver either. Here, it is done mainly for the sake of pedagogy to help students understand the interplay between a mathematical algorithm that looks good on paper and its hardware implementation which is not as promising as one would expect.

Jacobi and Gauss-Seidel wouldn't be solvers I'd even consider for a real problem. But they are simple enough that anyone with a basic understanding of linear algebra and programming can follow along. But much research requires me to run simulations on thousands (if not hundreds of thousands) of cores, and there, off-the-shelf solver implementations will often not cut it.