Poisson's Discrete Equation for 2 dimensions with Interfaces I'm trying to solve the Poisson's Discrete Equation when there are interfaces. I'm solving it for the electromagnetic potential for a given grid. When developing my equations, I've found the following relation:
\begin{align}\label{eq:EDP_x}
    \mathbf{A}_{i,j} = \frac{1}{2}
    \Bigg \{
    & \mathbf{A}_{i,j+1} \bigg[1 - \frac{\mu_{i, j+1}}{\mu_{i, j+1} + \mu_{i, j-1}} \bigg] +
    \mathbf{A}_{i,j-1} \frac{\mu_{i, j+1}}{\mu_{i, j+1} + \mu_{i, j-1}} \\ 
    +& \mathbf{A}_{i+1,j} \bigg[1  - \frac{\mu_{i+1, j}}{\mu_{i+1, j} + \mu_{i-1, j}} \bigg] +
    \mathbf{A}_{i-1,j} \frac{\mu_{i+1, j}}{\mu_{i+1, j} + \mu_{i-1, j}}  +
    \frac{\Delta x^{2}}{2} f_{i,j} \Bigg \}
\end{align}
Where $\mathbf{A}_{i,j}$ is the electromagnetic potential for a given $(i, j)$, $\mu_{i, j}$ is the permeability at that same point. Since the coefficients that multiply the neighbours of the target point vary for each point, I don't know exactly how to fit my equation to the known solution (the one displayed in the wikipedia's article). Is solving this equation matricially even possible? If so, how can I generate the left-side matrix, which is composed (in the general case) of the coefficients and the identity? Thank you very much in advance!
 A: Let's group the solution $\boldsymbol{A}$ into a vector (either row or column-major ordered). Let's do row-major here:
$$ \boldsymbol{A} = \begin{pmatrix} A_{1,1} \\ A_{1, 2} \\ \vdots  \\ A_{1, N_x} \\ A_{2, 1} \\ \vdots \\ A_{2, N_x} \\ \vdots \\ A_{N_y, N_x} \end{pmatrix}$$
Then, of course, you can set up a linear system of the form
$$
M (\mu) \boldsymbol{A} = \frac{\Delta x^2}{2} \boldsymbol{f}$$
to compute $\boldsymbol{A}$.
$$ M(\mu) = \begin{pmatrix} 2 & 1 - \frac{\mu_{i+1,j}}{\mu_{i+1,j} + \mu_{i-1,j} } & \dots & 1 - \frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1} } & & & \dots \\
\frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1}}& 2 & 1 - \frac{\mu_{i+1,j}}{\mu_{i+1,j} + \mu_{i-1,j} } & \dots & 1 - \frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1} } \\ 
& \ddots & \ddots & \dots & & \\
\frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1}} & \dots & \frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1}}& 2 & 1 - \frac{\mu_{i+1,j}}{\mu_{i+1,j} + \mu_{i-1,j} } & \dots & 1 - \frac{\mu_{i,j+1}}{\mu_{i,j+1} + \mu_{i,j-1} } \\
\ddots & & \ddots& & \ddots & & \ddots\end{pmatrix}$$
