Can I take advantage of a nearly banded A in AX=b?

Question

I am working on a 1D drift-driffusion problem in a finite-difference (FD) approach.

I hade 3 equations per node ($3N$ in total): electron continuity $E_i$, Poisson $P_i$, hole continuity $H_i$. With corresponding variables: electron concentration $n_j$, potential $V_j$, hole concentration $p_j$.

By rearranging the problem has a banded matrix, the Jacobian $A$ takes the following form, with shape $(3N\times 3N)$

$$A = \begin{pmatrix} A_{1,1} & A_{1,2} & & & & & \\ A_{2,1} & A_{2,2} & A_{2,3} & & & & \\ & A_{3,2} & A_{3,3} & A_{3,4} & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & A_{N-2,N-3} & A_{N-2,N-2} & A_{N-2,N-1} & \\ & & & & A_{N-1,N-2} & A_{N-1,N-1} & A_{N-1,N} \\ & & & & & A_{N,N-1} & A_{N,N} \\ \end{pmatrix}$$

where

$$ A_{i, j} = \begin{pmatrix} {\partial E_i}/{\partial n_j} & {\partial E_i}/{\partial V_j} & {\partial E_i}/{\partial p_j} \\ {\partial P_i}/{\partial n_j} & {\partial P_i}/{\partial V_j} & {\partial P_i}/{\partial p_j} \\ {\partial H_i}/{\partial n_j} & {\partial H_i}/{\partial V_j} & {\partial H_i}/{\partial p_j} \end{pmatrix} $$

Owing to that data structure, I have been using scipy.linalg.solve_banded to solve for $X$ in

$$AX = b$$

where $b$ is a $(3N \times 1)$ shaped vector.

My problem comes from that I want to add an equation that controls the potential at one end (say on the left) based on the value of total fluxes of $n$ and $p$ at both ends.

As a result the additional equation depends on $(n_1, V_1, p_1, n_2, V_2, p_2)$ but also $(n_{N-1}, V_{N-1}, p_{N-1}, n_N, V_N, p_N)$.

Therefore, if I was to include this new equation in $A$ to form a $(3N+1 \times 3N+1)$ Jacobian $A'$, the bandwith would suddenly become much larger. Though I can optimize this by inserting the new equation in the middle row of $A$, the banded representation would remain very sparse.

My question is then: is there a way I can solve for $X'$ in $$A' X' = b'$$

while keeping the advantages of the banded structure of $A$?

---

EDIT addressing comments

Here is a visual example of the Jacobian before (left) and after reordering (right). The last column/line, highlighted in yellow, is the additional equation that fails to fit the banded structure. Elements in dark gray background are non-zero, elements in white are zeros, elements in light gray background are zero, but a future version of my physical could in principle make them non-zero so an ideal solution would not exclude that possibility.

(click to enlarge)

Disclaimer: I do not pretend to be the author of this ordering method. I have read it a few years ago in a software manual I believe ...