For the Neumann problem $\,(\ast)\,$ in a bounded domain $U\subset\mathbb{R}^n$, $n\geqslant 2$, satisfying the cone condition, to prove that assumption
$$
f\in \{ L^2(U)\,\colon\;\int\limits_{U}f\,dx=0\}\tag{1}
$$
implies the existence of a weak solution $u\in H^1(U)$, it is convenient to introduce the space
$$
\widetilde{H}^1(U)=\{w\in H^1(U)\colon\,\int\limits_{U}\!w\,dx=0\}.
$$
Notice that $\widetilde{H}^1(U)$ is a Hilbert space with inner product
$$
(u,v)\overset{\rm def}{=}\int\limits_{U}\nabla u\cdot\nabla v\,dx
$$
satisfying the condition
$$
(u,u)=0\;\;\Longrightarrow\;\;u=0
$$
by virtue of the Poincaré inequality
$$
\|u\|^2_{L^2(U)}\leqslant C\int\limits_{U}|\nabla u|^2\,dx
\quad \forall\,u\in \widetilde{H}^1(U)\tag{2}
$$
which requires certain regularity of the boundary $\partial U$. Note that the cone condition is not
precisely the regularity of $\partial U$ for $(2)$ to be valid — it just proves to be the least complicated suitable general restriction on $\partial U$. Denote
$$
\bar{u}\overset{\rm def}{=}\frac{1}{|U|}\int\limits_{U}u\,dx,
$$
with notation $|U|$ standing for the $n$-dimensional Lebesgue measure of domain
$U\subset \mathbb{R}^n$. Since $u-\bar{u}\in \widetilde{H}^1(U)$ for any
$u\in H^1(U)$, the Poincaré inequality can be as well rewritten in the form
$$
\|u-\bar{u}\|^2_{L^2(U)}\leqslant C\int\limits_{U}|\nabla (u-\bar{u})|^2\,dx
=C\int\limits_{U}|\nabla u|^2\,dx
\quad \forall\,u\in H^1(U).
$$
The rest of the proof is easy. Consider a linear functional
$$
\Lambda(v)=\int\limits_{U}fv\,dx
$$
on $\widetilde{H}^1(U)$. Due to $(2)$, the linear functional $\Lambda$ is bounded
on the Hilbert space $\widetilde{H}^1(U)$. Hence, by the Riesz representation theorem, there is a unique $u\in\widetilde{H}^1(U)$ such that
$$
\Lambda(v)=(u,v)\quad \forall\,v\in \widetilde{H}^1(U),\tag{3}
$$
which immediately implies the integral identity
$$
\int\limits_{U}\nabla u\cdot\nabla v\,dx=\int\limits_{U}fv\,dx
\quad \forall\,v\in \widetilde{H}^1(U).\tag{4}
$$
To complete the proof, notice that, in fact, $(4)$ is valid as well for all
$u\in H^1(U)$. Indeed, due to the assumption $(1)$, for any $v\in H^1(U)$ we have
$$
\int\limits_{U}fv\,dx=\int\limits_{U}f(v-\bar{v})\,dx=
\int\limits_{U}\nabla u\cdot\nabla (v-\bar{v})\,dx=
\int\limits_{U}\nabla u\cdot\nabla v\,dx
$$
by virtue of $(3)$ since $v-\bar{v}\in \widetilde{H}^1(U)$. Thus, there is a unique
$u\in\widetilde{H}^1(U)\subset H^1(U)$ such that
$$
\int\limits_{U}\nabla u\cdot\nabla v\,dx=\int\limits_{U}fv\,dx
\quad \forall\,v\in H^1(U).
$$
Q.E.D
Remark. Being valid for general real bilinear forms, not necessarily symmetric, the Lax-Milgram theorem looks too much advanced for this rather trivial case when all the inner product axioms are met by the symmetric bilinear form
$\,(\cdot,\cdot)$. Generally, the Lax-Milgram theorem is to be applied in cases where
the Riesz representation theorem is inapplicable, e.g., in case of a Dirichlet problem for the equation $-\Delta u+\partial_{x_m}u=f$.
