Solve this specific large sparse system of linear equations
The key is the fact that $C_i = - I$, i.e., the appropriately size identity matrix is trivially nonsingular and can be used to nullify the blocks above each copy. We have the following row equivalence $$ A = \begin{pmatrix} B_0 & B_1 & & & \\ & -I & B_2 & & \\ & & -I & B_3 & \\ & & & -I & B_4 \\ D & & & & -I \\ \end{pmatrix} \sim \begin{pmatrix} B_0 & B_1 & & & \\ & -I & B_2 & & \\ & & -I & B_3 & \\ B_4D & & & -I & \\ D & & & & -I \\ \end{pmatrix} \sim \begin{pmatrix} B_0 & B_1 & & & \\ & -I & B_2 & & \\ B_3B_4D & & -I & & \\ B_4D & & & -I & \\ D & & & & -I \\ \end{pmatrix}\sim \begin{pmatrix} B_0 & B_1 & & & \\ B_2B_3B_4D & -I & & & \\ B_3B_4D & & -I & & \\ B_4D & & & -I & \\ D & & & & -I \\ \end{pmatrix}\sim \begin{pmatrix} B_0+B_1B_2B_3B_4D & & & & \\ B_2B_3B_4D & -I & & & \\ B_3B_4D & & -I & & \\ B_4D & & & -I & \\ D & & & & -I \\ \end{pmatrix} $$ You can now factor the block $B_0 + B_1B_2B_3B_4D$ and compute $x_0$. The rest will be just be forward substitution, but compute them in descending order $x_4$, $x_3$, $x_2$, $x_1$ so that you do not have to regenerate the spike on the left.
You could try formulating it as a feasibility optimization problem such as $$\begin{equation} \begin{array}{rrclcl} \displaystyle \min_{x} &0 \\ \textrm{s.t.} & A x & = & b \\ \end{array} \end{equation}$$
and feed it into a linear programming solver based on a branch-and-price algorithm such as GCG. Such algorithms are usually designed to exploit structured variables via a Dantzig-Wolfe decomposition. In fact, GCG is a branch-price-and-cut solver, meaning it also adds cutting planes in addition to the column generation approach via DW-decomposition. Very often such structured problems like yours can be solved fast via such a decomposition algorithm.
In the above problem formulation you can use any objective function as you are only looking for a feasible solution of your system $Ax=b$. Technically, you could replace the 0-objective function with any other objective function but why not just keep it simple?
The easiest quite fast way is probably to use Carl-Christians algebraic solution.
If you are interested in low-level calculational tricks for these matrices you can keep reading below line.
Given the matrix structure there are some key matrix elements at the tips of the triangles. If you let the solution diffuse from there as in forward-substitution-like scheme or implicitly with CG you will gain a lot.
Maybe you can get some inspiration from another question of mine, if I can find it. Any triangular matrix you can factor into sparse di-diagonal matrix "in-place" multiplication. This saves calculations. But I don't remember where I wrote this question right now so I will try explain.
You can factor: $$B_i = F_1F_2\cdots F_m$$
Where each $F_i$ is an identity matrix except for the entries (i,i) and (i,i+1) which need to be determined.
Whooops $m$ matrices you think, but $m$ is huge, but since we only need store 2 values for each, it will only need $2m$ in total for each $B_i$, compared to $m^2/2$ which the tridiagonal $B_i$ ones need.
In practice for $m=10^5$ you could reduce memory fingerprint from tens gigabytes to singles of megabytes.