Given the complexity of analytical solutions, it is helpful to numerically approximate solutions to steady state, constant heat conduction, two dimensional conduction.

$$ \frac{\partial^2T}{\partial x^2}+\frac{\partial^2T}{\partial y^2}=0 $$

One numerical method of solving our two dimensional conduction equation is the finite difference method. This method is very similar to our generalized technique for solving linear systems of temperatures, but is able to handle heat generation through the material.

The first considerations are general system parameters, prescriptions for boundary conditions, and meshing settings.

# === General Parameters ===
w_x = 1.                            # m, Width in the x-direction
w_y = 1.                            # m, Width in the y-direction
k = 1.                              # W/m/k, Conductivity function
gd = 0.1                            # W/m^3, Heat generation

# === Boundary Conditions ===
T_N = lambda x : 100*4*x*(w_x-x) + 100  # C, Boundary temperature function
Qd_S = lambda x : 0                     # W, Boundary insulation
Qd_E = lambda y : 50*y*(w_x-y)          # W, Boundary temperature gradient function
T_W = lambda y : 100                    # C, Boundary temperature

# === Meshing Settings ===
delx = 0.05             # m, Node spacing in the x-direction
dely = delx             # m, Node spacing in the y-direction
nx = int(w_x/delx+1)    # Number of columns
ny = int(w_y/dely+1)    # Number of rows
x = lambda ix : ix*delx # Function for resolving x-position from x-direction index
y = lambda iy : iy*dely # Function for resolving y-position from y-direction index

We will construct a linear system out of a governing equation for each node, for which we will initialize two matrices, with which we will solve for our temperatures. In order to build this system, we must iterate over every node in the mesh. We will traverse along each row in the mesh, starting from the top left, moving right across the first row, and then repeating for each successive row.

# === System Initialization and Nodal Traversal ===
T_sys_A = np.zeros((nx*ny, nx*ny))
T_sys_B = np.zeros((nx*ny,))
for iy in range(ny):        # Traversal across rows
    for ix in range(nx):    # Traversal across columns
        row = iy
        col = ix
        num = iy * nx + ix          # Index of the node at row iy and column ix
        num_N = (iy - 1) * nx + ix  # Index above the current node
        num_S = (iy + 1) * nx + ix  # Index below the current node
        num_E = iy * nx + (ix + 1)  # Index right of the current node
        num_W= iy * nx + (ix - 1)   # Index left of the current node

Within each step of this traversal, our governing of equation must be configured differently based on the whether the node is an interior node or various types of boundary conditions. Note that the order of your if statements will change which boundary conditions govern the corners.

# === Node Type Evaluation ===
        if(row == 0):               # North boundary condition
            # ...
        elif(row == ny - 1):        # South boundary condition
            if(col == 0):           # South-east corner condition
            # ...
            elif(col == nx - 1):    # South-west corner condition
            # ...
            else:                   # South edge condition
            # ...
        elif(col == nx - 1):        # East boundary condition
            # ...
        elif(col == 0):             # West boundary condition
            # ...
        else:                       # Interior node
            # ...

At each node, based on the identified node type, we must configure our system matrices to express the appropriate governing equation. The most general form of conservation of energy is as follows, stating that the conduction in to a node plus the heat generation in the node is equivalent to the negative of any external heat transfer. To derive all the following equations, this is expanded by the definition of conduction and then simplified. The following derivations assume $\Delta x=\Delta y$, that $k,\dot g$ are invariant, but in principle the technique works even without these assumptions.

$$ \dot Q_N+\dot Q_S+\dot Q_E+\dot Q_W+\dot gV_i=-\dot Q_{b~i} $$

$$ k_(\Delta x\frac{T_N-T}{\Delta y}+ \Delta x\frac{T_S-T}{\Delta y}+ \Delta y\frac{T_E-T}{\Delta x}+ \Delta y\frac{T_W-T}{\Delta x})+\Delta x\Delta y\dot g=-\dot Q_{b~i} $$

Specified Value

$$ (1)T_i=T_b $$

Heat Flux from the North (Edge Node)

$$ (4)T_i+(-2)T_S+(-1)T_E+(-1)T_W=(\frac{2\dot Q_b+\Delta x\Delta y\dot g}{k}) $$

Heat Flux from North and West (Corner Node)

$$ (4)T_i+(-2)T_{Si}+(-2)T_{Ei}=(\frac{4\dot Q_b+\Delta x\Delta y\dot g}{k}) $$

Interior Node

$$ (4)T_i+(-1)T_N+(-1)T_S+(-1)T_E+(-1)T_W=(\frac{\Delta x\Delta y\dot g}{k}) $$