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.