How to can I transform the 2D cuasi Laplace equation with variable coefficients to finite difference scheme? I want to solve $$\frac{\partial}{\partial x}\left(\frac{1}{\rho(x,y)}\frac{\partial \Phi}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{1}{\rho(x,y)}\frac{\partial \Phi}{\partial y}\right)= 0$$ in a rectangular domain and with Neumann boundary conditions. I want to applied $\frac{\delta^2_x \phi}{\delta_x \rho}+ \frac{\delta^2_y \phi}{\delta_x \rho}=0$, with $\delta_x$ and $\delta^2_x$ operators of finite differences schemes. I mean, I'm not sure of using the $\delta$ operator up and down to discretize the Laplace equation with variable terms.
$$\delta^2_x \phi = \frac{\phi_{i+1,j}-2\phi{i,j}+\phi_{i+1,j}}{\Delta x^2}$$
$$\delta_x \phi = \frac{\phi_{i+1,j}-\phi_{i-1,j}}{2\Delta x}$$
Suggestions are welcome.
 A: If I'm not mistaken, you want do
$$
\frac{\partial}{\partial x} \left(\frac{1}{\rho} \frac{\partial \phi}{\partial x}\right) \simeq 
\frac{\delta_x^2 \phi}{\delta_x \rho} = 
\frac{\frac{\phi_{i+1, j} - 2 \phi_{i+1, j} + \phi_{i-1, j}}{\Delta x^2}}{\frac{\rho_{i+1, j} - \rho_{i-1, j}}{2 \Delta x}}.
$$
I think this is not a good idea. It you observe,
$$
\frac{\delta_x^2 \phi}{\delta_x \rho} \simeq \frac{\partial_{xx}\phi}{\partial_x \rho},
$$
and you have
$$
\partial_{x}\left(\frac{1}{\rho}\partial_x \phi\right) = \frac{\rho \, \partial_{xx} \phi - \partial_x \rho \, \partial_x \phi}{\rho^2},
$$
if $\rho(x,y)$ is differentiable. Then you could discretize the above form. However, there is another way to approximate the operators more directly (and I'm think it is better). We need to define the finite difference operator
$$
\delta^*_x f(x_i, y_j) = \frac{f(x_i+\Delta x/2 \,,\, y_j) - f(x_i-\Delta x/2 \,,\, y_j)}{\Delta x} = \frac{f_{i+\frac{1}{2} \,,\, j} - f_{i-\frac{1}{2} \,,\, j}}{\Delta x}.
$$
We approximate each differential operator $\partial_x$ by the finite difference operator $\delta^*_x$.
Let be $\alpha(x, y) = 1/\rho$, we have
$$
\partial_{x}\left(\alpha(x_i, y_j) \partial_x \phi(x_i, y_j)\right) \simeq \delta^*_x \left(\alpha_{i,j} \, \delta^*_x \phi_{i,j}\right)
= \frac{\alpha_{i+1/2, \,j} \, \delta^*_x \phi_{i+1/2,\,j} -  \alpha_{i-1/2, \,j} \, \delta^*_x \phi_{i-1/2,\,j} }{\Delta x} = \frac{\alpha_{i+1/2, \,j} \, (\phi_{i+1,\,j}-\phi_{i,j}) -  \alpha_{i-1/2, \,j} \, (\phi_{i,j}-\phi_{i-1,\,j})}{\Delta x^2} 
= \frac{\alpha_{i+1/2, \,j} \, \delta^*_x \phi_{i+1/2,\,j} -  \alpha_{i-1/2, \,j} \, \delta^*_x \phi_{i-1/2,\,j} }{\Delta x} = \frac{\alpha_{i+1/2, \,j} \, \phi_{i+1,\,j} - (\alpha_{i+1/2, \,j} + \alpha_{i-1/2, \,j}) \phi_{i,j} +  \alpha_{i-1/2, \,j} \, \phi_{i-1,\,j}}{\Delta x^2}.
$$
To approximate $\alpha_{i+1/2, \,j}$, you can take the arithmetic mean,
$$
\alpha_{i+1/2, \,j} = \frac{\alpha_{i+1, \,j}+\alpha_{i, \,j}}{2}.
$$
