I am considering a system that can be simplified to the following problem. Let $\Omega$ be a bounded open subset of $\mathbb{R}^n$ with smooth boundary, $n \geq 3$, and consider the following coupled system of PDEs for a pair of positive functions $u, v$ on $\Omega$ : $$ \left\lbrace \begin{array}{rl} -\Delta u + u &= |\nabla v|^2 u^{-N-1},\\ -\Delta v + v &= u^N \end{array} \right. $$ together with the Neumann boundary conditions $$ \left\lbrace \begin{array}{rl} \partial_\nu u &= 0,\\ \partial_\nu v &= u^N \end{array} \right. $$ Here $N = \frac{2n}{n-2}$ but the exact value is probably not relevant. I can prove the existence of a weak solution to this problem, namely $u$ is such that $u^{N/2+1} \in W^{1, 2}(\Omega, \mathbb{R})$ (multiply the equation for $u$ by $u^{N+1}$ and integrate by parts over $\Omega$) but I am struggling to improve the regularity of the solution...
What can be done, by means of the sub- and super solutions, is to prove that, if $\nabla v \in L^q(\Omega, \mathbb{R})$ for some $q > 1$ then $u \in W^{2, q}(\Omega, \mathbb{R})$. The difficulty is that we need Sobolev regularity for $u^N$ near the boundary to prove regularity for $v$ and raising $u$ to the power $N$ causes a drastic loss of regularity...
Any help would be invaluable.
For those who are interested in the original problem, please see https://arxiv.org/abs/1403.5655 (Section 4).
