Double beta decay (DBD) is a rare nuclear process of great interest due to its potential to provide information about physics beyond the Standard Model (BSM). For example, the discovery of the neutrinoless double-beta (0 νββ ) decay mode could give answers to fundamental issues about lepton number conservation, CP and Lorentz symmetry violation in the weak sector, and about still unknown neutrino properties such as: are neutrinos Dirac or Majorana particles?; neutrino absolute masses; what is the correct hierarchy of the neutrino masses?; are there sterile neutrinos?, etc. Theoretically, the first important stage in the DBD study consists in the derivation of half-lives formulas and precisely computation of the nuclear matrix elements (NME) and phase space factors (PSF) entering these formulas, for different decay modes and transitions to ground or excited states of the daughter nuclei. Reliable computations of these quantities result in reliable predictions of DBD half-lives and constrains of the BSM parameters related to the possible mechanisms that may contribute to the 0 νββ decay. In this paper we give first a short review of the theoretical challenges in the computation of the NME and PSF for double-beta decay. Then we present a new, more consistent, approach to calculate these quantities, namely to calculate directly their product.