Your discretization of $\nabla \cdot (k \nabla T)$ is incorrect. An immediate correction is to replace $\Delta x$ and $\Delta y$ in the 2nd & 5th lines with $(\Delta x)^2$ and $(\Delta y)^2$ -- in other words,
\begin{align*}
(\nabla \cdot (k \nabla T))_{i,j} &= k(x_i, y_j) \frac{T_{i+1,j} + T_{i-1,j} - 2T_{i,j}}{(\Delta x)^2} + k_x(x_i, y_j) \frac{T_{i+1,j} - T_{i,j}}{\Delta x} \\
& \qquad \cdots + k_y(x_i, y_j) \frac{T_{i,j+1} - T_{i,j}}{\Delta y} + k(x_i, y_j) \frac{T_{i,j+1} + T_{i,j-1} - 2T_{i,j}}{(\Delta y)^2}. \tag{*}
\end{align*}
Note, however, that $\frac{T_{i+1,j} + T_{i-1,j} - 2T_{i,j}}{(\Delta x)^2}$ is a second-order accurate approximation of $T_{xx}$, while $\frac{T_{i+1,j} - T_{i,j}}{\Delta x}$ is only a first-order accurate approximation of $T_x$. Therefore, a numerical scheme based on (*) can be expected to have at most first-order accuracy (in the $L^\infty$ norm) in general. For a second-order accurate discretization, see ($\dagger$) at the end. In order to motivate it, however, let's focus on 1D, where $\nabla \cdot (k \nabla T) = (k T_x)_x$.
A second-order discretization of $(k T_x)_x$ is
\begin{align*}
((k T_x)_x)_i &= \frac{k(x_{i+1/2}) (u_x)_{i+1/2} - k(x_{i-1/2}) (u_x)_{i-1/2}}{2(h/2)} \tag{1} \\
&= \frac{k(x_{i+1/2}) \frac{u_{i+1} - u_i}{2(h/2)} - k(x_{i-1/2}) \frac{u_i - u_{i-1}}{h}}{2(h/2)}, \tag{2}
\end{align*}
which reduces to
\begin{align*}
(\nabla \cdot (k \nabla T))_i &= \frac{k(x_{i+1/2}) u_{i+1} + k(x_{i-1/2}) u_{i-1} - (k(x_{i+1/2}) + k(x_{i-1/2})) u_i}{h^2}, \tag{3}
\end{align*}
with $x_{i+1/2} = \frac12 (x_i + x_{i+1})$.
How did this calculation work?
- The right-hand side of (1) is simply the central finite difference (FD) discretization of $\nu_x$ where $\nu = kT_x$. In contrast with the forward- or backward- FD discretizations of $\nu_x$, which are first-order accurate, a central FD discretization of $\nu_x$ is second-order accurate.
- Equation (2) is a repeat -- it's the central FD discretization of $u_x$.
- Intuitively, since we've used second-order accurate approximation, we might claim that (3) is second-order accurate -- one can prove that this is true, for example, by computing the local truncation error.
- Note: The denominators of (1) and (2) involve "$2(h/2)$" -- rather than $h$ -- since the denominator (of the central FD discretization of $\nu_x$ or $u_x$) is $= 2\times\text{step size}$, where -- in this case -- $\text{step size} = h/2$.
Equation (3) is easily extended to 2D in a direction-by-direction manner, e.g.,
\begin{align*}
(\nabla \cdot (k \nabla T))_{i,j} &= \frac{k(x_{i+1/2,j}) u_{i+1,j} + k(x_{i-1/2,j}) u_{i-1,j} - (k(x_{i+1/2,j}) + k(x_{i-1/2,j})) u_{i,j}}{h^2} \\
& \qquad \cdots + \frac{k(x_{i,j+1/2}) u_{i,j+1} + k(x_{i,j-1/2}) u_{i,j-1} - (k(x_{i,j+1/2}) + k(x_{i,j-1/2})) u_{i,j}}{h^2}, \tag{$\dagger$}
\end{align*}
with $x_{i+1/2,j}$, $x_{i,j+1/2}$, etc., defined analogously.
