I'm sorry I ignore the condition of $a_{ij}$, but there is a Moser's iteration for $-\Delta u + V(x)u=0$. For instance, see Struwe's Variational Methods Appendix B lemma B.3. And I think this iteration doesn't need the continuous of $a$.
More precisely, integral the equation $Lu+cu=0$ by $u\min\{|u|^{2s}, L^2\}\in H_{0}^{1}(\Omega)$, then we have $$\int a_{ij}u_i u_j \min\{|u|^{2s}, L^2\} +\frac{s}{2}\int _{|u|^s\leq L} a_{ij} (|u|^2)_i (|u|^2)_j|u|^{2s-2}\leq c\int |u|^2\min\{|u|^{2s}, L^2\}.$$ For $s$ small enough such that $u\in L^{2+2s}(\Omega)$ we have $$\int a_{ij}u_i u_j \min\{|u|^{2s}, L^2\} \text{ is bounded uniformly with} ~L.$$ Hence let $L \to \infty$, we have $$\int |D|u|^{s+1}|^2<\infty,$$ which means, by Sobolev embedding, $u\in L^{\frac{(2s+2)n}{n-2}}(\Omega)$. By iteration, you can obtain that $u\in L^{q}(\Omega)$ for all $q>1$.