Consider the fallowing result:
Assume $a^{ij} \in C^{1}(U)$, $f \in L^{2}$ and let $u \in H^{1}(U)$ be a weak solution of $$ -\sum_{i,j=1}^{n}(a^{ij}u_{x_i})x_j = f, \text{ in }U $$ Then $u \in H^{2}_{loc}(U)$.
I would like to know if there's any result like this with $a^{ij} \in L^{\infty}(U)$. Any reference with regularity results with coefficients in $L^{\infty}(U)$ would be very helpful!
Yes results similar to this are proved in Elliptic Partial Differential Equations of Second Order by Gilbarg and Trudinger. Chapter 8 is devoted to the study of such equation, in particular see Section 8.9. The following is Theorem 8.24 in Section 8.9 which I've adjusted to your situation.
Suppose that $a^{ij}$ are measurable on a domain $\Omega$ and there exists positive constants $\lambda,\Lambda$ such that $$a^{ij}(x) \zeta_i \zeta_j \geqslant \lambda \vert \zeta \vert^2 \qquad \text{for all }x\in \Omega, \, \zeta \in \mathbb R^n $$ and $$\sum_{i,j} \vert a^{ij}(x) \vert^2 \leqslant \Lambda^2 \qquad \text{for all }x\in \Omega. $$
You may also find the discussion in Section 2.5 of Regularity Theory for Elliptic PDE by Fernadez-Real and Ros-Oton interesting. Here they consider the case $a^{ij}$ are continuous. In particular, they mention that you cannot expect weak solutions to your PDE to have a bounded $C^1$ norm - only $C^{1-\varepsilon}$ for all $\varepsilon>0$.