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?

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).