Based on Legendre wavelets, a new efficient numerical integration method is proposed to estimate double integrals. This new method is expressed in terms of the operational integration matrices. Those which were introduced by Parsian, are calculated in an approximate way. In this work, we also propose the exact computation of these matrices. Better precision of the new method is observed during a comparative numerical study with the composite Simpson method.