Poisson equation in non-homogenous media If the dielectric constant be different in x, y, and z directions, how the Poisson equation would change. Is it right to say:
$$ \epsilon_{xx} d^2V/dx^2 + \epsilon_{yy} d^2V/dy^2  + \epsilon_{zz} d^2V/dz^2  = -\rho$$
 A: Linear anisotropic materials can be modeled as having a field $P^i = \epsilon_0 ~ \chi^i_j ~ E^j$. This follows two conventions: the first is "Einstein summation rules"; you sum over any index which is repeated both upper and lower. The second convention is that upper indices correspond to a vector space while lower indices correspond to its dual space. This difference is huge when you're talking about crystals, where the atoms have a certain set of real-space vectors $\hat b_i$ which tell you how the crystal is laid out, so any vector is $\vec v = \hat b_i v^i$, while the dual-space vectors perpendicular to them $\hat b^i$ are defined so that $\hat b^i \cdot \hat b_j = \delta^i_j$ is the Kronecker delta (1 if $i = j$, 0 otherwise). If you work through it you will find that there is also a matrix called the metric $g^{ij}$ which relates the two components via $v^i = g^{ij} v_j.$
So we see that we replace the normal electric susceptibility $\chi$ with a susceptibility tensor $\chi^i_j$. This causes $D^i = \epsilon^i_j ~E^j$ for $\epsilon^i_j = \epsilon_0 (\delta^i_j + \chi^i_j)$. This means that if $E_k = -\partial_k \phi$ for some potential $\phi$, the relevant equation for a bulk with no free charges is $$\partial_i D^i = \partial_i (\epsilon^i_j ~ g^{jk} ~ \partial_k\phi) = (\partial_i \epsilon^{ik}) \partial_k \phi + \epsilon^{ik} \partial_i \partial_k \phi = 0.$$ 
For a homogeneous ($\chi$ not depending on position) linear anisotropic substance, we know $\partial_n \epsilon^{ik} = 0$ and so we can confidently discard that first term.
The second term is trickier, but there will be certainly some basis which diagonalizes $\epsilon^i_j = e^{i} ~ \delta^i_j$, and it will usually be the crystal basis, so you will still have the expression $\epsilon^j ~ g^{jk} ~ \partial_j ~ \partial_k \phi ~=~ \epsilon^j ~ \partial_j ~ \partial^j \phi = 0.$ In general it looks like this expression will "tweak" the $\epsilon^j$ coefficients with the metric before they enter into an equation substantially similar to yours.
If the material is non-homogeneous then of course the first term in the above expression is some vector field $t^k = \partial_i (\epsilon^i_j g^{jk})$ which also contributes a $(\vec t \cdot \nabla) \phi$ term to the mix.
