A new variational method for studying the equilibrium states of an interacting particles system has been proposed. The statistical description of the system is realized by means of a density matrix. This method is used for description of conduction electrons in metals. An integral equation for the electron distribution function over wave vectors has been obtained. The solutions of this equation have been found for those cases where the single-particle Hamiltonian and the electron interaction Hamiltonian can be approximated by a quite simple expression. It is shown that the distribution function at temperatures below the critical value possesses previously unknown features which allow to explain the superconductivity of metals and presence of a gap in the energy spectrum of superconducting electrons.