We present and examine an approximate but ab initio many-body approach, viz., potential harmonics expansion method (PHEM), which includes two-body correlations for dilute Bose-Einstein condensates. Comparing the total ground state energy for three trapped interacting bosons calculated in PHEM with the exact energy, the new method is shown to be very good in the low density limit which is necessary for achieving Bose-Einstein condensation experimentally.