Suppose we have a smooth solution to a system of PDE (for $\vec{u}$) given by $L\vec{u} = \vec{f}$ on the unit ball (no boundary conditions given). Assume that the coefficients and $\vec{f}$ are smooth functions and that this system is elliptic (i.e., the symbol of L is invertible and is uniformly so, in the sense the its determinant is bounded away from 0).


I have $C^0$ bounds on the coefficients of $L$ and $C^{k}$ (for any k) bounds on $\vec{f}$. What interior a priori estimates can one get on $\vec{u}$ ? 



A priori estimates for elliptic systems with bounded coefficients
In general, estimates for elliptic systems with bounded coefficients are only slightly better than the energy class that solutions live in. Consider the case 
$$L({\bf u}) = \text{div}(A(x)D{\bf u}) = \partial\_i(A^{ij}\_{\alpha\beta}(x) u^{\beta}\_j) = 0.$$ 
Here ${\bf u} = (u^1,\,...,\,u^m)$ is a map from $\mathbb{R}^n \rightarrow \mathbb{R}^m$, and $\lambda |p|^2 \leq A^{ij}\_{\alpha\beta}(x)p^{\alpha}\_ip^{\beta}\_j < \Lambda |p|^2$ for all $p \in M^{m \times n}$ and all $x$. Such ${\bf u}$ are local minimizers of the energy $\int A(D{\bf u},\,D{\bf u})\,dx$, so $W^{1,\,2}$ estimates are expected. By using energy minimality at all scales one can get the improvement
$$\left(\int\_{B\_{1/2}} |D{\bf u}|^{2 + \delta}\,dx \right)^{\frac{2}{2+\delta}} < C\int\_{B\_1} |{\bf u}|^2\,dx$$
for some $\delta > 0$ small depending on the ellipticity constants. This is known as the ''reverse-Holder theory'' (see e.g. the book of Giaquinta, ''Multiple Integrals in the Calculus of Variations and Nonlinear Elliptic Systems.") 


In the scalar case $m = 1$ solutions are $C^{\alpha}$ by the De Giorgi-Nash theorem. The key is the maximum principle; in the scalar case it is never energetically favorable to go outside the boundary values, so the oscillation of the solution decreases as we focus near a point.


In the vectorial case, the above energy estimate is optimal. Heuristically, it is not always energetically favorable to truncate a component, so a map might pay less by going outside the boundary data.


To show optimality one can play with modifications of the well-known De Giorgi example (also in Giaquinta's book). To see the idea, let $n = m$ and consider the discontinuous map $\nu = \frac{x}{|x|}$. Then $D\nu = r^{-1}P$ where $P$ is the matrix that projects tangent to to the sphere. Since $P$ is perpendicular to the matrix $N$ that projects in the radial direction, $\nu$ minimizes the energy with coefficients $N \otimes N$. Since these coefficients are degenerate elliptic, we perturb the them slightly:
$$A(x) = \delta I\_{n^2} + (N + \epsilon rD\nu) \otimes (N + \epsilon rD\nu).$$
Imposing $\text{div}(AD\nu) = 0$ gives (after some calculation) 
$$(n-1)\delta = \epsilon(n-1)(n-2-\epsilon(n-1)),$$
so the coefficients can be made uniformly elliptic provided $n > 2$ and $\epsilon$ small.


Using coefficients of the same form, one can show that $|x|^{-\alpha}\nu$ are homogeneous solutions to uniformly elliptic systems in $\mathbb{R}^n \backslash \{0\}$ for all $\alpha \neq \frac{n-2}{2}$. The homogeneity $-\frac{n-2}{2}$ is exactly $W^{1,\,2}$ critical; for $\alpha$ smaller than this value, these maps are in $W^{1,\,2}$, and one can in fact construct smooth approximations to the coefficients and the map that are solutions. Furthermore, as $\alpha \rightarrow \frac{n-2}{2}$ we see that the energy estimate is optimal.


(Interestingly, the only solutions homogeneous of degree $-\frac{n-2}{2}$ 
in $B\_1 \backslash \{0\}$ are trivial. Indeed, multiplying the equation by ${\bf u}$ and integrating by parts gives
$$\int\_{B\_1 \backslash B\_{\epsilon}} |D{\bf u}|^2\,dx < C\int\_{\partial(B\_1 \backslash B\_{\epsilon})} A(D{\bf u}, {\bf u}\otimes \nu)\,ds.$$ 
Since $D{\bf u}\cdot {\bf u}$ is $1-n$-homogeneous, the right side is bounded independent of $\epsilon$, and by $W^{1,\,2}$ criticality, the left side blows up as $\epsilon \rightarrow 0$ unless ${\bf u}$ is constant.)


