The question is basically in the title. I want to evolve an image $u_0$ (i.e. a $n_1\times n_2$ resolution set of discrete values in $[0,1)$, which is a discrete function on $\{0,\ldots,n_1-1\}\times\{0,\ldots,n_2-1\}$) by the PDE $$\partial_tu=\nabla\cdot\kappa(t,\nabla u)\nabla u\tag1$$ with Neumann boundary conditions. $\kappa$ is real-valued. I want to simulate $(1)$ using the Euler scheme together with a spatial discretization of the space derivatives.
Now, my problem is that there are various possibility how the spatial derivatives could be approximated by forward, backward or central finite differences. For the boundary condition as well. Since we are working with discrete images, I really wonder which of these I should use and which spacing (probably $1$?) I should use.
I mostly worry about stability. I think a bad scheme might explode, even though it should not. Any help is highly appreciated!