Problem statement: Due to the long range nature of interactions of the N-body systems, direct computation of the Coulomb potential energy involves O(N2) operations. To decrease such complexity, a simple Multilevel Summation method has been developed. Approach: In the frame of the Multilevel Summation method, the two-body interaction is decomposed into two parts: a local part and a smooth part. The local part vanishes beyond some cut-off distance; hence, its contribution to the potential energy is calculated in O(N) operations. In contrast to some common fast summation methods, the smooth part is calculated in real space on a sequence of grids with increasing meshsize in O(N) operations. Results: The method is tested on the calculation of the Madelung constants of ionic crystals in one, two and three dimensional cases. For a cut-off distance equals three times the meshsize of the ionic crystal, an error less than 0.01% is obtained. Conclusion: In computing the coulomb lattice sums of charge systems consisting of N bodies, the Multilevel Summation method decreases the complexity to O(N) operations