How to discretize the bioheat equation to solve with finite difference method. 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}
 A: 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$.
