An extended Debye-Hiickel theory with fourth order gradient term is developed for electrolyte solutions; namely, the electric potential phi(r) of the bulk electrolyte solution can be described by del(2)phi(r) = kappa(2)phi(r) + L-Q(2)del(4)phi(r), where the parameters kappa and L-Q 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.