I need to discretize the bioheat equation in 2D to solve it numerically, and I'm choosing the finite difference method to solve it.
In steady state situation
$$\rho_t c_t \frac{\partial T}{\partial t} = 0$$
But I don't know how to discretize
$ \overbrace{k_t \left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right)}^{\substack{\text{heat conduction term } }} - \underbrace{\rho_b c_b w_b (T_{(x,y)}-T_b) }_{\substack{\text{advection by perfusing} \\ \text{blood (cooling source)}}} + \overbrace{Q_m}^{\substack{\text{metabolic} \\ \text{heating}}} + \underbrace{P}_{\substack{\text{energy dissipation} }} = 0 $
For a simple heat conduction problem in steady state the discretization can be as:
$$ T_{i,j} = \frac{1}{4}\left[T_{i+1,j} + T_{i-1,j} + T_{i,j+1} + T_{i,j-1}\right] $$
but this equation has 3 other terms, that I don't know how to include into the discretization.
UPDATE I have made this discretization, but I'm not sure if i's well done just discretizing the derivative term, and isolating the $T_{i,j}$.
The discretization over 2D grid system of the heat conduction term using finite difference method, is:
\begin{equation} k \left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right) = k \left[ \frac{T_{i-1,j} - 2T_{i,j} + T_{i+1,j}}{\Delta x^2} + \frac{T_{i,j+1} - 2T_{i,j} + T_{i,j-1}}{\Delta y^2} \right] \end{equation}
so:
\begin{equation} k \left[ \frac{T_{i-1,j} - 2T_{i,j} + T_{i+1,j}}{\Delta x^2} + \frac{T_{i,j+1} - 2T_{i,j} + T_{i,j-1}}{\Delta y^2} \right] - \rho_b c_b w_b (T_{i,j}-T_b) + Q_m + P_m = 0 \end{equation}
Using an equally spaced grid system in $x = y$ axis:
\begin{equation} d = \Delta x = \Delta y \end{equation}
Then the term $T_{i,j}$ has to be isolated in order to get the expresion of the temperature in the node ${i,j}$.
\begin{equation} k \left[ \frac{T_{i-1,j} - 2T_{i,j} + T_{i+1,j} + T_{i,j+1} - 2T_{i,j} + T_{i,j-1}}{d^2} \right] - \rho_b c_b w_b (T_{i,j}-T_b) + Q_m + P_m = 0 \end{equation}
\begin{equation} k \left[ \frac{T_{i-1,j} + T_{i+1,j} + T_{i,j+1} + T_{i,j-1} - 4T_{i,j}}{d^2} \right] - \rho_b c_b w_b (T_{i,j}-T_b) + Q_m + P_m = 0 \end{equation}
\begin{equation} k \left[ \frac{T_{i-1,j} + T_{i+1,j} + T_{i,j+1} + T_{i,j-1} - 4T_{i,j}}{d^2} \right] = \rho_b c_b w_b (T_{i,j}-T_b) - Q_m - P_m \end{equation}
\begin{equation} \frac{T_{i-1,j} + T_{i+1,j} + T_{i,j+1} + T_{i,j-1} - 4T_{i,j}}{d^2} = \frac{\rho_b c_b w_b (T_{i,j}-T_b) - Q_m - P_m}{k} \end{equation}
\begin{equation} T_{i-1,j} + T_{i+1,j} + T_{i,j+1} + T_{i,j-1} - 4T_{i,j} = d^2 \frac{\rho_b c_b w_b (T_{i,j}-T_b) - Q_m - P_m}{k} \end{equation}
\begin{equation} \boxed{ T_{disc} = T_{i-1,j} + T_{i+1,j} + T_{i,j+1} + T_{i,j-1} }\end{equation}
\begin{equation} T_{disc} - 4T_{i,j} = d^2\left[\frac{\rho_b c_b w_b (T_{i,j}-T_b) - Q_m - P_m}{k}\right] \end{equation}
\begin{equation} - 4T_{i,j} = d^2 \left[ \frac{\rho_b c_b w_b (T_{i,j}-T_b) - Q_m - P_m}{k} \right] - T_{disc} \end{equation}
\begin{equation} - 4T_{i,j} = \frac{d^2 \rho_b c_b w_b (T_{i,j}-T_b) - d^2 Q_m - d^2 P_m }{k} - T_{disc} \end{equation}
\begin{equation} - 4T_{i,j} = \frac{d^2 \rho_b c_b w_b (T_{i,j}-T_b) - d^2 Q_m - d^2 P_m - kT_{disc}}{k} \end{equation}
\begin{equation} - 4kT_{i,j} = d^2 \rho_b c_b w_b (T_{i,j}-T_b) - d^2 Q_m - d^2 P_m - kT_{disc} \end{equation}
\begin{equation} \boxed{ \alpha = \rho_b c_b w_b }\end{equation}
\begin{equation} - 4kT_{i,j} = d^2 \alpha(T_{i,j}-T_b) - d^2 Q_m - d^2 P_m - kT_{disc} \end{equation}
\begin{equation} - 4kT_{i,j} = d^2 \alpha T_{i,j} - d^2 \alpha T_b - d^2 Q_m - d^2 P_m - kT_{disc} \end{equation}
\begin{equation} - d^2 \alpha T_{i,j} - 4kT_{i,j} = - d^2 \alpha T_b - d^2 Q_m - d^2 P_m - kT_{disc} \end{equation}
\begin{equation} d^2 \alpha T_{i,j} + 4kT_{i,j} = d^2 \alpha T_b + d^2 Q_m + d^2 P_m + kT_{disc} \end{equation}
\begin{equation} T_{i,j}(d^2 \alpha + 4k) = d^2 \alpha T_b + d^2 Q_m + d^2 P_m + kT_{disc} \end{equation}
\begin{equation} T_{i,j}(d^2 \alpha + 4k) = d^2(\alpha T_b + Q_m + P_m) + kT_{disc} \end{equation}
\begin{equation} \boxed{ T_{i,j} = \frac{d^2(\alpha T_b + Q_m + P_m) + kT_{disc}}{(d^2 \alpha + 4k)}} \end{equation}
You correctly discretized the Laplace operator with a five-point stencil but I'm not sure what you are trying to achieve by solving for $T_{ij}$. There are different ways to go about but the standard approach is the following.
Assume you have $N_x, N_y$ cells in $x$ and $y$ direction, respectively. You then construct a one-dimensional vector
$$ T = \begin{pmatrix} T^{[1]} \\ \vdots \\ T^{[N_y]} \end{pmatrix} $$
where $T^{[i]} = \begin{pmatrix} T_{i,1} & \cdots & T_{i,N_x} \end{pmatrix}^T$ represents the $i$-th row of your grid.You then collect all terms containing the various $T_{ij}$ etc on the left-hand side and all other (source) terms on the right-hand side. You then end up with a matrix equation
$$ A T = f $$
where $A$ has the form
$$ A =\begin{pmatrix} V & I & & & & & \\ I & V & I & & & \\ & I & V & I & & \\ & & \ddots & \ddots & \ddots & & \\ & & & I & V & I \end{pmatrix} $$
with $V$ being a tri-diagonal matrix (representing $\partial_x^2$) and $I$ being a diagonal matrix (representing the $y$ derivatives). The diagonal in $V$ consists of all $T_{ij}$ in your equation.
Likewise, all terms which do not contain any $T_{k,l}$ are collected in $f$. You will also need to incorporate the boundary conditions in $f$. This is easy for Dirichlet boundary conditions.
You then typically use an iterative matrix solver to find your solution vector $T$.