Matrix has icomplete rank, so the system is solvable only if
$$
\operatorname{rank} M = \operatorname{rank} \begin{pmatrix}M \;\big|\; f\end{pmatrix}.
$$
Due to truncation error, the last may not hold (but would hold, if you're using a conservative approximation, I suppose).
You can use QR decomposition to deal with that. Suppose $M = QR$ where $Q$ is orthogonal and $R$ is upper triangular matrix. The rank of $R$ is equal to the rank of $M$ and that is one less than number of its rows. So
$$
R = \begin{pmatrix}
\star & \star & \dots & \star & \star\\
0 & \star & \dots & \star & \star\\
0 & 0 & \ddots & \star & \star\\
0 & 0 & \dots & \star & \star\\
0 & 0 & \dots & 0 & 0\\
\end{pmatrix} = \begin{pmatrix}R' & h\\0&0\end{pmatrix}
$$
Solving the system $Mx = f$ is equivalent to solving the system
$$
Rx=Q^\top f \equiv \begin{pmatrix}z'\\\alpha\end{pmatrix}, \alpha \in \mathbb R.
$$
For the system to have a solution the last element $\alpha$ of $Q^\top f$ should be zero. Even if it's not, let's fix it as zero. Let's also fix the last element of $x$ as zero, since it can be any number.
$$
x = \begin{pmatrix}x'\\0\end{pmatrix}
$$
Now we have a nonsingular triangular system
$$
R' x' = z'
$$
which can easily be solved.
Note, that for the scheme you're using matrix $M$ is sparse, so you need sparse version of $QR$ decomposition. Also, since $Q$ is not sparse for QR even if $M$ was, you need to perform QR for extended matrix $(M|f)$ to compute both $R = Q^\top M$ and $Q^\top f$ at the same time.
Appendix
I've implemented this method in Matlab and also the method of one equation elimination. For correct problems they give the same solution.
For testing purposes I used square domain $[-1,1]\times[-1,1]$ and the following problem
$$
\Delta u = 0\\
\frac{\partial u}{\partial n} = g(x, y), \quad (x,y) \in \partial G.
$$
This problem has well known solvability criterion, that is
$$
0 = \left(\int_G \Delta u dx dy = \int_{\partial G} \frac{\partial u}{\partial n} d\Gamma\right) = \int_{\partial G} g d\Gamma
$$
Due to numerical errors it never can be exactly achieved, so from exact mathematical point of view, the numerical problem is always ill-posed. So we have to deal with regularized problem, that is to solve a problem with
$$
\int_{\partial G} g d \Gamma \neq 0.
$$
That's where the methods differ. They eliminate the inconsistency in various ways.
Equation elimination
When we are using equation elimination, that is removing one governing equation and replacing it with
$$
u_{i,j} = 0
$$
for example, we state that all the inconsistency should eliminated via $(i,j)$ point. In fact, the real problem we are solving now is
$$
\Delta u = \delta(x - x_i, y - y_i) \int_{\partial G} g d\Gamma\\
\frac{\partial u}{\partial n} = g(x, y), \quad (x,y) \in \partial G.
$$
so we put a source or a sink at $x_i, y_i$ to resolve the inconsistency. For the test problem I want $u(x,y) = \frac{x^2}{2} - \frac{y^2}{2}$ to be the solution, so
$$
g(-1, y) = 1,\quad g(1, y) = 1\\
g(x, -1) = -1, \quad g(x, 1) = -1
$$
Do demonstrate the effect of inconsistency I would take $g(x, -1) = 0.5$.
Note the pole-like thing at the middle - that was the point with $u_{i,j} = 0$ equation. At every other point the function satisfies $\Delta u = 0$ and the boundary conditions.
QR method
Again, let
$$
Mu = f
$$
be decomposed as
$$
Mu \equiv QRP^\top u = f\\
RP^\top u = Q^\top f.
$$
I've added permutation matrix $P$ to allow column reordering of the matrix $M$ which improves sparse QR a lot.
Let again,
$$
R = \begin{pmatrix}
R' & h\\
0 & 0
\end{pmatrix}, \quad Q^\top f = \begin{pmatrix}
z'\\\alpha
\end{pmatrix}
$$
So we state that
$$
x = \begin{pmatrix}
R^{-1} z'\\
0
\end{pmatrix}
$$
is a regularized solution to
$$
Rx = z \Leftrightarrow Mx = f
$$
since $Mx = f$ simply does not have a solution. But what is the real problem we are solving now?
Recall that $M$ is singular, we also know its left zeroing vector. Let $\omega = \frac{1}{\|\omega\|}(1, 1, \dots, 1)^\top$. Then
$$
\omega^\top M
$$
is simply the sum of every row of $M$ and that is zero (actually, it is implementation-dependent, but the idea remains the same). So the last row of $R$ is nothing more than $\omega^\top M$, and $\alpha = \omega^\top f \sim \int_{\partial G} g d\Gamma$ (up to some constant multiplication).
If we plug the regularized solution $x$ to the system $Rx = z$ we observe a residual:
$$
R x \neq z,\quad R x = z + \begin{pmatrix}0\\-\alpha\end{pmatrix}
$$
and by multiplying by $Q$ we have
$$
Mx = f - \alpha \omega
$$
since $\omega$ is exactly the last column of $Q$. So now we are solving a problem where inconsistency is smoothed over the whole domaing $G$ and its border $\partial G$, roughly speaking
$$
\Delta u = -\frac{1}{|G|} \int_{\partial G} g d\Gamma
$$
This method applied to the same problem gives the following solution
It looks fine, but, actually, the equation is violated (slightly) at every node and boundary conditions are also violated (also slightly).
Finally, there's no much difference what method to use if you're solving a well-posed problem, like arising from
$$
\nabla \phi = \mathbf E \Rightarrow \Delta \phi = \operatorname{div} \mathbf E.
$$
The Matlab code I've used you can find here
Best Answer
The meaning of $H(u)$ is a matter of convention. A Neumann boundary condition (BC) writes $\frac{\partial u}{\partial n}(x,y) = k(x,y)$. The derivative with respect to $n$ means "in the direction normal to the boundary" with an agreed convention on the sign of this normal. I'll get back to this later. For now, you can put $k$ on the left-hand side and call that $H(u)$, i.e. $H(u)=\frac{\partial u}{\partial n}(x,y) - k(x,y)$ . In your example I suspect $k\equiv 0$, but let's assume it is not for now.
The Neumann BC involves a derivative and we need to represent it using finite differences. You should use an order of finite difference that is the same as the one you are using inside of the domain. The 5 points stencil is second order. We thus need to express the Neumann BC using a second order scheme. We have two main options: central differences or forward/backward differences.
To go further we need to know which boundary we are talking about. It seem your example concern a boundary in a plane $x=cst$ which is the "inlet" of the domain. Such a boundary would typically correspond to the value $j=0$. Note that all the things derived from now on are specific to this boundary. The normal to this boundary, pointing towards the inside of the domain, is thus $n=e_y$. The derivative "in the direction normal to the boundary", will be taken as the derivative along $y$ that is we need to ensure that $\frac{\partial u}{\partial y}(x_i,y_j) = k(x_i,y_j)$, where $j=0$. Now we need to represent that using a second order finite difference scheme: central differences or forward differences (forward since we are at the inlet and we can't really use information backward). Let's take the example of the central differences.
If we use central differences, we will go backward by one point. I've said we couldn't use information backward, but in this case it introduces an "unknown" that we call a ghost point, and we can solve for its value since we know $k$. The central difference approximation writes: \begin{align} j=0, \quad \frac{\partial u}{\partial y}(x_i,y_j)=k(x_i,y_j) \rightarrow \frac{U_i^{j-1}-U_i^{j+1}}{2 \Delta y} = k(x_i,y_j) \end{align} which can be solved for the ``ghost point'' as \begin{align} j=0,\quad U_i^{j-1} = U_i^{j+1} + 2 k(x_i,y_j)\Delta y \end{align} In your case $\Delta y=h$. Note that if we were at the outlet (say $j=n$), the ghost point would be at $n+1$. At this point, the normal towards the domain is along $-e_y$ and we will have: \begin{align} j=n, \quad U_i^{j+1} =U_i^{j-1} + 2 k(x_i,y_j)\Delta y \end{align} This is a matter of convention between what we put in the variable $k$ and what we choose to be the normal.
Now we need to ensure that the boundary condition is met for the Poisson equation. We write the Poisson equation at the boundary point itself (that's just the general formula, at $j=0$): $$ j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+U_{i}^{j+1}+[U_{i}^{j-1}]-4U_{i}^{j}) = f_i^j $$ we can replace $[U_i^{j-1}]$ in this formula by the value we obtained above: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+U_{i}^{j+1}+[ U_i^{j+1} + 2 k(x_i,y_j)h]-4U_{i}^{j}) = f_i^j$$ Gathering the terms: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+2U_{i}^{j+1}+ 2 k(x_i,y_j)h-4U_{i}^{j}) = f_i^j$$ You see that the term $U_{i}^{j+1}$ appears twice indeed. Now in your case, I suspect $k=0$ so the equation at the boundary is: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+2U_{i}^{j+1}-4U_{i}^{j}) = f_i^j$$
The formula above is close to the formula you give. Yet, there is no reason for $U_{i+1}^{j}+U_{i-1}^{j}$ to be equal to zero. The only case where this can happen is if these two points are on a boundary where a Dirichlet value of 0 is specified. In this case your domain is just made of three nodes in the $x$ direction, counting the boundary nodes.