Previous comments gave you good suggestions, I try to add some more.

Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.

If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include the values at boundaries to unknowns (you might done that already). For four nodes in the corners of your rectangle, you have to change in the corresponding row of matrix two -1 to -2, for all other rows corresponding to boundary nodes change only one -1 to -2, it must be exactly the single entry that is inside of rectangle so not on boundary. As noted by someone else previously, the sum off all entries for each row in the matrix must be equal to zero, that is another check for you if you are doing it correctly.

The idea behind finite difference discretization for zero Neumann boundary conditions is that you imagine that you have one row or one column of nodes next to each side of your rectangle with additional (artificial) unknowns. It means that you may use the standard stencil with 4 on the diagonal and four time -1 offdiagoals. Additional unknowns can be eliminated by adding discrete equations that correspond to zero Neumann boundary condition, e.g. \partial_x \varphi = 0, that takes in the discrete form (hopefully the notation is clear) e.g. $(\varphi_{1 j} - \varphi_{-1 j})/dx = 0$, so the artificial unknowns are equal $\varphi_{-1 j} = \varphi_{1 j}$. Plug this into the discrete equation for unknowns $\varphi_{0 j}$ and you get what I did suggest. By the way this corresponds to second order accurate approximation that you should prefer.

Another point is that applying zero Neumann boundary condition for Poisson equations gives you underdetermined system where there are formally infinitely many solutions, but they differ only by constant (it means if you subtract two solution functions you get constant function). Once you fix one value e.g. to 0 (as proposed by someone else), you choose one of those infinitely many solution, so your problem has then unique solution. It is like you define Dirichlet boundary condition (e.g. equals 0) in one single point at boundary. Of course, then such node must be excluded from your unknowns, so you do not have a row in your matrix for this node.

Sometimes your linear solver can take care for this, but this is another story, the suggestion above may be much simpler.

Previous comments gave you good suggestions, I try to add some more.

Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.

If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include the values at boundaries to unknowns (you might done that already). For four nodes in the corners of your rectangle, you have to change in the corresponding row of matrix two -1 to -2, for all other rows corresponding to boundary nodes change only one -1 to -2, it must be exactly the single entry that is inside of rectangle so not on boundary. As noted by someone else previously, the sum off all entries for each row in the matrix must be equal to zero, that is another check for you if you are doing it correctly.

The idea behind finite difference discretization for zero Neumann boundary conditions is that you imagine that you have one row or one column of nodes next to each side of your rectangle with additional (artificial) unknowns. It means that you may use the standard stencil with 4 on the diagonal and four time -1 offdiagoals. Additional unknowns can be eliminated by adding discrete equations that correspond to zero Neumann boundary condition, e.g. $\partial_x \varphi = 0$, that takes in the discrete form (hopefully the notation is clear) e.g. $(\varphi_{1 j} - \varphi_{-1 j})/dx = 0$, so the artificial unknowns are equal $\varphi_{-1 j} = \varphi_{1 j}$. Plug this into the discrete equation for unknowns $\varphi_{0 j}$ and you get what I did suggest. By the way this corresponds to second order accurate approximation that you should prefer.

Another point is that applying zero Neumann boundary condition for Poisson equations gives you underdetermined system where there are formally infinitely many solutions, but they differ only by constant (it means if you subtract two solution functions you get constant function). Once you fix one value e.g. to 0 (as proposed by someone else), you choose one of those infinitely many solution, so your problem has then unique solution. It is like you define Dirichlet boundary condition (e.g. $\varphi =0$) in one single point at boundary. Of course, then such node must be excluded from your unknowns, so you do not have a row in your matrix for this node.

Sometimes your linear solver can take care for this, but this is another story, the suggestion above may be much simpler.

Previous comments gave you good suggestions, I try to add some more.

Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.

If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include the values at boundaries to unknowns (you might done that already). For four nodes in the corners of your rectangle, you have to change in the corresponding row of matrix two -1 to -2, for all other rows corresponding to boundary nodes change only one -1 to -2, it must be exactly the single entry that is inside of rectangle so not on boundary. As noted by someone else previously, the sum off all entries for each row in the matrix must be equal to zero, that is another check for you if you are doing it correctly.

The idea behind finite difference discretization for zero Neumann boundary conditions is that you imagine that you have one row or one column of nodes next to each side of your rectangle with additional (artificial) unknowns. It means that you may use the standard stencil with 4 on the diagonal and four time -1 offdiagoals. Additional unknowns can be eliminated by adding discrete equations that correspond to zero Neumann boundary condition, e.g. \partial_x \phi = 0, that takes in the discrete form (hopefully the notation is clear) e.g. $(\phi_{1 j} - \phi_{-1 j})/dx = 0$, so the artificial unknowns are equal $\phi_{-1 j} = \phi_{1 j}$. Plug this into the discrete equation for unknowns $\phi_{0 j}$ and you get what I did suggest. By the way this corresponds to second order accurate approximation that you should prefer.

Another point is that applying zero Neumann boundary condition for Poisson equations gives you underdetermined system where there are formally infinitely many solutions, but they differ only by constant (it means if you subtract two solution functions you get constant function). Once you fix one value e.g. to 0 (as proposed by someone else), you choose one of those infinitely many solution, so your problem has then unique solution. It is like you define Dirichlet boundary condition (e.g. equals 0) in one single point at boundary. Of course, then such node must be excluded from your unknowns, so you do not have a row in your matrix for this node.

Sometimes your linear solver can take care for this, but this is another story, the suggestion above may be much simpler.

Previous comments gave you good suggestions, I try to add some more.

Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.

If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include the values at boundaries to unknowns (you might done that already). For four nodes in the corners of your rectangle, you have to change in the corresponding row of matrix two -1 to -2, for all other rows corresponding to boundary nodes change only one -1 to -2, it must be exactly the single entry that is inside of rectangle so not on boundary. As noted by someone else previously, the sum off all entries for each row in the matrix must be equal to zero, that is another check for you if you are doing it correctly.

The idea behind finite difference discretization for zero Neumann boundary conditions is that you imagine that you have one row or one column of nodes next to each side of your rectangle with additional (artificial) unknowns. It means that you may use the standard stencil with 4 on the diagonal and four time -1 offdiagoals. Additional unknowns can be eliminated by adding discrete equations that correspond to zero Neumann boundary condition, e.g. \partial_x \varphi = 0, that takes in the discrete form (hopefully the notation is clear) e.g. $(\varphi_{1 j} - \varphi_{-1 j})/dx = 0$, so the artificial unknowns are equal $\varphi_{-1 j} = \varphi_{1 j}$. Plug this into the discrete equation for unknowns $\varphi_{0 j}$ and you get what I did suggest. By the way this corresponds to second order accurate approximation that you should prefer.

Another point is that applying zero Neumann boundary condition for Poisson equations gives you underdetermined system where there are formally infinitely many solutions, but they differ only by constant (it means if you subtract two solution functions you get constant function). Once you fix one value e.g. to 0 (as proposed by someone else), you choose one of those infinitely many solution, so your problem has then unique solution. It is like you define Dirichlet boundary condition (e.g. equals 0) in one single point at boundary. Of course, then such node must be excluded from your unknowns, so you do not have a row in your matrix for this node.

Sometimes your linear solver can take care for this, but this is another story, the suggestion above may be much simpler.