In general, the estimates for systems are just 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 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 
$$(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.)