For the Laplace Equation with fixed values of $\phi$ at the boundary, under a simple electrostatic design with an homogeneous material, $$ \nabla^2\phi=0\\ \phi(S)=\phi_0(S), \tfrac{\partial\phi}{\partial n}(S)=\phi'_0(S) $$ the following finite differences iterative step worked fine for me, $$ \hat\phi(x,y)={1 \over 4} (\phi(x-1,y)+\phi(x+1,y)+\phi(x,y-1)+\phi(x,y+1)) $$ including this additional step for efficiency (ref.): $$ \phi(x,y)=k(\hat\phi(x,y)-\phi(x,y))+\phi(x,y)\\ k=2-2\pi/n \\ $$
The Faraday and Gauss Maxwell' Equations turn into the Laplace Equation if $\epsilon=\epsilon_0$, $$ \nabla\times E=0 \to -\nabla\times \nabla\phi=0 \\ \nabla\cdot\epsilon E=0 \to -\nabla\cdot\epsilon \nabla \phi=0 $$
but if I have a variable material $\epsilon=\epsilon(x,y)$, (which under almost no cases should consist in discrete steps or others discontinuities, but only in smooth variations) I do not longer have the Laplace Equation, and the previous rule is no longer valid, since the equation in this case is: $$ \nabla\cdot\epsilon\nabla\phi=0 $$
A first trial to satisfy this new equation, is following the same finite differences and try to obtain a new iteration step, $$ \nabla\cdot\epsilon\nabla\phi=\nabla\epsilon\cdot\nabla\phi+\epsilon\nabla^2\phi=0 $$ so following the same approach with differences, a possible iterative step $\hat{\hat\phi}$ could perhaps be: $$ \chi(x+1,y)=\epsilon(x+1,y)/\epsilon(x,y)-1\\ \hat{\hat\phi}(x,y)={\chi(x+1,y)\phi(x+1,y)+\chi(x,y+1)\phi(x,y+1)+4\hat\phi(x,y) \over 4+\chi(x+1,y)+\chi(x,y+1)} $$ but I did not find any proper reference to validate it, for going further with this.
This problem is not at all new, but I have not found a consistent reference to check, as I do not know how to deal with this, instead making trial and error while people perhaps already solved it the proper and smart way.
Question is, how should I solve the Laplace Equation for a Non Homogeneous Material.
