We have that $$\tag{1}\int\nabla u_\epsilon\nabla\varphi+\int\beta(u_\epsilon)\varphi=f\varphi,\ \forall\ \varphi\in H_0^1$$
By taking $\varphi=u_\epsilon$, we conclude $$\tag{2}\int|\nabla u_\epsilon|^2\leq\int fu_\epsilon$$
which implies by Holder inequality that $\|\nabla u_\epsilon\|_2$ is bounded, therefore, we can suppose that $u_\epsilon \rightharpoonup u$ in $H_0^1$. Moreover, we can suppose that $u_\epsilon \to u$ in $L^2$. We combine these convergences with $(2)$ to get $$\tag{3}\int |\nabla u|^2\leq \int fu$$
On the other hand, define $\mathcal{K}=\{w\in H_0^1:\ w\geq 0\}$. As you can easily verify, $\int\beta(u_\epsilon)w\leq 0$ in $\mathcal{K}$, hence, from $(1)$, we conclude that $$\tag{4}\int\nabla u_\epsilon\nabla\varphi\geq \int f\varphi,\ \forall\ \varphi\in \mathcal{K}$$
Again, we use the convergences we have, to conclude from $(4)$ $$\tag{5}\int\nabla u\nabla\varphi\geq\int f\varphi,\ \forall\ \varphi\in\mathcal{K}$$
We conclude from $(3)$ and $(5)$ the desired inequality. To conclude that $u\in\mathcal{K}$ and $u$ is unique, you can argue like this:
I - $u$ is the minimizer of some functional $I:\mathcal{K}\to\mathbb{R}$, what is $I$?
II - In this particular setting, for $u$ to minimizes $I$ is equivalently for $u$ to satisfies the variational inequality.
III - The solution of the optimization problems is unique.
As @RayYang suggested this book is a good one to understand it better, however I would like to point out that the main argument here can be found in any good book of convex analysis.
It is worth to note that what we are proving here is that $-I'(u)\in \mathcal{N}_{\mathcal{K}}(u)$, where $\mathcal{N}_{\mathcal{K}}(u)$ is the normal cone of $u$ with respect to $\mathcal{K}$ and when $I$ is a convex differentialbe function, II is valid, i.e. $-I'(u)\in \mathcal{N}_{\mathcal{K}}(u)$ if and only if $u$ minimizes $I$.
Best Answer
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$.