I have question about the discretization of the poisson equation $$\triangle u =f$$ on a two-dimensional grid as described in https://en.wikipedia.org/wiki/Discrete_Poisson_equation.
In this example, they use the natural ordering, in other words $u_{1,1},u_{2,1},...$ etc and create an equation system of the form: $$A\cdot U=b$$ where $U=(u_{1,1},u_{2,1},\dots,u_{1,2},\dots)^{T}$
but remove the boundarie values $u_{0,0},u_{1,0},u_{0,1},...$etc.
This works fine if I have dirichlet boundary conditions like $u(x,y)=g(x,y)$ for some function on the boundary, but what if I have Neumann-Condition, so that the boundary values depend on discritized points $u_{i,j}$ which needs to be calculated too? One approach was to modify $A$, by modifying entries of the matrix $A$ so that I don't need to add $u_{0,1},u_{1,0},\dots$ etc to the vector $U$. Another approach was to increase the dimension of the matrix $A$ and add the boundary values in the vector $U$. BUt by doing so, I would loose properties like symmetry of Matrix A.
My question is if anyone knows a paper how one would implement non-zero boundary conditions?