We study the regularity of weak solutions to the elliptic system in divergence form divA(x, Du)=0 in an open set Ω of R^n, n≥2. The vector field A(x.ξ), A: Ω×R^(m×n)→R^(m×n), has a variational nature in the sense that A(x,ξ)= Dξ f(x,ξ), where f:Ω×R^(m×n)→R is a convex Carathéodory integrand; i.e., f=f (x,ξ) is measurable with respect to x∈R^n and it is a convex function with respect to ξ∈R^(m×n). If m=1 then the system reduces to a partial differential equation. In the context m>1 of general vector-valued maps and systems, a classical assumption finalized to the everywhere regularity of the weak solutions is a modulus-dependence in the energy integrand; i.e., we require that f(x,ξ)=g(x,|ξ|), where g:Ω×[0,∞)→[0,∞) is measurable with respect to x∈ R^n and it is a convex and increasing function with respect to the gradient variable t∈[0,∞).