I am not sure where's the best place to ask this question but I am wondering if anyone has some insight on my problem.
I am trying to numerically solve the following PDE in divergence form:
$\nabla \cdot(c \nabla u) = 0$ subject to the boundary conditions:
$x \in [0,1] \times[0,2]\subset \mathbb{R}^2$, $u(0,x_2) = 0, u(1,x_2) = 1$ and $u_{x_2}(x_1,0) = u_{x_2}(x_1,2) = 0$.
I'm trying out very simple functions $c(x_1,x_2)$.
Suppose that he coefficient has the form $c(x_1,x_2) = 2 + (8x_1^2 - 8x_1 + 1)(2x_2 - 1) +(2x_1 - 1)(8x_2^2 - 8x_2 + 1)$.
The graph of such a coefficient is as follows:

However, when I use Matlab's FEM solver even with the finest mesh, the solution looks like this:

I find it so strange that such a smooth coefficient $c(x_1,x_2)$ can yield such irregular of a solution.
But when I apply scaling parameters to my field $c(x_1,x_2)$, i.e. it now becomes
$c(x_1,x_2) = 2 + 0.2(8x_1^2 - 8x_1 + 1)(2x_2 - 1) + 0.1(2x_1 - 1)(8x_2^2 - 8x_2 + 1)$, the solution is well-behaved and respects the boundary conditions:
Does anyone have an idea of what's going on? I do not have the intuition behind this...

Your diffusion coefficient $c$ is negative at places in your domain. For example $c(0,2)=-12$. The PDE is ill-posed when $c$ changes sign, since the equation is then not elliptic.
Your scaled diffusion coefficient is positive everywhere on your domain, so all is good.