Let $K$ be a compact set in $\mathbb R^n,$ $f\in C^\infty(K),$ and $c:\mathbb R \to \mathbb R$ be a smooth function. Consider an element $u\in H_0^1(K),$ satisfying the PDE $ -\Delta u + c(u) = f$ in weak sense, i.e.
$$ \int_K \nabla u\cdot \nabla v+\int_K vc(u)=\int_K fv. $$
I think we can establish something similar to $u \in C^\infty(K).$ Indeed, we can consider a standard elliptic regularity argument: let $v=-D^{-h}_kD^{h}_ku$, where $D^{h}_ku(x)=\frac{u(x+he_k)-u(x)}{h}$. Then, we can estimate: $$ \int \nabla u\cdot \nabla v=\|D^{h}_k \nabla u\|_{L^2(K)}^2. $$ The difficult part is the estimate on the nonlinearity. Firstly I try to get $u\in H^2.$ We can easily write $$ \int_K (f-c(u))D^{-h}_kD^{h}_ku\leq (\|f\|_{L^2} + \|c\circ u\|_{L^2}) \|D^{-h}_kD^{h}_k u\|_{L^2} $$ If we can show that $\|c\circ u\|_{L^2}<\infty,$ then it is just a matter of simple calculations to confirm the gain in regularity. However, this is difficult. In dimensions $n=1,2$, we have the continuous Sobolev embedding $H^1 \to L^\infty$, so we can treat $u$ as a bounded function and hence $\|c\circ u \|_{L^\infty}<\infty,$ and since $K$ has finite measure, $c\circ u\in L^2$. So the desired regularity is obtained. However, in higher dimensions, we do not have such nice embedding.
Do we have (at least) $H^2$ regularity in higher dimensions?
Apparently it is more realistic to have a bound on the growth of $c$. However, I find out that this question is asked in Lawrence C. Evan's Partial Differential Equations book on page 366 without any bounds on $c$.