In a bounded Lipschitz domain in R n , we consider a second-order strongly elliptic system with symmetric principal part written in divergent form. We study the Neumann boundary value problem in a generalized variational (or weak) setting using the Lebesgue spaces H σ p (Ω) for solutions, where p can differ from 2 and σ can differ from 1. Using the tools of interpolation theory, we generalize the known theorem on the regularity of solutions, in which p = 2 and |σ − 1| < 1/2, and the corresponding theorem on the unique solvability of the problem (Savaré, 1998) to p close to 2. We compare this approach with the nonvariational approach accepted in numerous papers of the modern theory of boundary value problems in Lipschitz domains. We discuss the regularity of eigenfunctions of the Dirichlet, Neumann, and Poincaré-Steklov spectral problems.