Dirichlet Boundary Condition finite difference method using sparse-matrix $Ax = b$ system

I am trying to solve the boundary value problem for heat equation:

$$
u_{xx} + u_{yy} = f(x,y)
$$

where the solution $u(x,y) in [0,1] times [0,1]$ and the Dirichlet boundary condition $u(x,y) = u_0$.

So I discretized the problem with $N$ interior points with step size $h = 1/(N+1)$

$$
w_{i,j-1} + w_{i+1, j} + w_{i-1,j} – 4 w_{i,j} + w_{i, j+1} = h^{2} f(x_{i}, y_{j})
$$

where $i,j in {1,2,…,N}$ and the boundary conditions are:

$$
w_{0,j} = w_{N+1, j} = u_{0} quadquadquad w_{i,0} = w_{i,N+1} = u_{0}
$$

where $i,j$ is from $1$ to $N$

I reduce this problem into sparse-matrix system. I know how to make the sparse matrix $A$ for this problem by using Python. The only issue is the RHS which is matrix $b$. For the case of $u_0 = 0$, I trivially evaluate the grid points at function $f(x,y)$. For the case $u_0 neq 0$, I do not know to retrieve matrix $b$. I tried to write down with the case $N = 3$ to see the general pattern as below (I put in column-order):

I hope anyone could help me understand how to retrieve that matrix $b$ in general case for Dirichlet boundary condition??. The code I wrote for this problem is in Python, but I want to understand the generality first.