An extended Debye–Hückel theory with fourth order gradient term is developed for electrolyte solutions; namely, the electric potential φ(r) of the bulk electrolyte solution can be described by ∇2φ(r) = κ2φ(r) + LQ2∇4φ(r), where the parameters κ and LQ are chosen to reproduce the first two roots of the dielectric response function of the bulk solution. Three boundary conditions for solving the electric potential problem are proposed based upon the continuity conditions of involving functions at the dielectric boundary, with which a boundary element method for the electric potential of a solute with a general geometrical shape and charge distribution is derived. Solutions for the electric potential of a spherical ion and a diatomic molecule are found and used to calculate their electrostatic solvation energies. The validity of the theory is successfully demonstrated when applied to binary as well as multicomponent primitive models of electrolyte solutions.