The longitudinal charge density of an electron beam in its equilibrium state is given by the solution of the Haissinski equation, which provides a stationary solution of the Vlasov-Fokker-Planck equation. The physical input is the longitudinal wake potential. We formulate the Haissinski equation as a nonlinear integral equation with the normalization integral stated as a functional of the solution. This equation can be solved in a simple way by the matrix version of Newtonss iteration, beginning with the Gaussian as a first guess. We illustrate for several quasi-realistic wake potentials. Convergence is extremely robust, even at currents much higher than nominal for the storage rings considered. The method overcomes limitations of earlier procedures, and provides the convenience of automatic normalization of the solution.