Discretize derivative of heat flux with variable conductivity How do i discretize this Temperature flow equation with the Fnite Difference method, when the conductivity K, is not constant?:
$$
    \frac{\partial}{\partial x}\left(K\frac{\partial T}{\partial x}\right) + 
    \frac{\partial}{\partial y}\left(K\frac{\partial T}{\partial y}\right) = 0  
$$
My attempt at this problem was to do the following:
$$
K(x)\frac{\partial^2 T}{\partial x^2} + 
    K(y)\frac{\partial^2 T}{\partial y^2} = 0 
$$
then discretize the second order partial equations:
$$
K(x) \left(\frac{T_{i+1,j} -2T_{i,j} + T_{i-1,j} }{\Delta x}\right) +        K(y) \left(\frac{T_{i,j+1} -2T_{i,j} + T_{i,j-1} }{\Delta y}\right) = 0
$$
However, i have been informed that this is an incorrect method as this working assumes a constant conductivity, where that is not the case in the problem. 
would someone be able to help me with this?
 A: Let us consider the 1D case for simplicity. The second-order centered finite-difference may be viewed as the composition of two first-order centered finite-differences:
$$
\left.\frac{\text d^2 T}{\text d x^2} \right|_{x=x_i}\simeq \frac{T_{i+1}-2T_{i}+T_{i-1}}{{\Delta x}^2} = \frac{1}{\Delta x}\left(\frac{T_{i+1}-T_{i}}{{\Delta x}} - \frac{T_{i}-T_{i-1}}{{\Delta x}}\right) .
$$
In the spatially-varying conductivity case, the latter finite-difference writes as
$$
\left.\frac{\text d}{\text d x}\left(K\frac{\text d T}{\text d x}\right)  \right|_{x=x_i}\simeq K_{i+1/2}\frac{T_{i+1}-T_{i}}{{\Delta x}^2} - K_{i-1/2}\frac{T_{i}-T_{i-1}}{{\Delta x}^2}\, ,
$$
where $K_{i\pm 1/2} = K\big(x_i\pm\frac{1}{2}\Delta x\big)$. It is incorrect to take $K$ out of the derivative as proposed in OP, since by the product rule, $(K T')' = K T'' + K'T'$ where prime $'$ denotes differentiation w.r.t. $x$.
This expression suggests to compute the finite-difference
\begin{aligned}
\left.\frac{\text d}{\text d x}\left(K\frac{\text d T}{\text d x}\right)  \right|_{x=x_i} &\simeq \frac{K_{i+1/2}+K_{i-1/2}}{2}\frac{T_{i+1}-2T_{i}+T_{i-1}}{{\Delta x}^2}\\
&\phantom{ \simeq } + \frac{K_{i+1/2}-K_{i-1/2}}{{\Delta x}}\frac{T_{i+1}-T_{i-1}}{2\,{\Delta x}}\,
\end{aligned}
which leads to the same method as before. Other suggestions and the same one can be found here.
In 2D, we have $\partial_x(K \partial_x T) + \partial_y(K \partial_y T)$ with conductivity $K(x,y)$ instead, but the same principle applies (see e.g. this post).
