A new model, based on the BCS approach, is especially designed to describe nuclear phenomena (A, Z ) → (A, Z ± 2) of double-charge exchange (DCE). Although it was proposed and applied in the particle-hole limit, by one of the authors [Krmpotić, Fizika B 14, 139 (2005)], it has not yet been applied within the BCS mean-field framework, nor has its ability to describe DCE processes been thoroughly explored. It is a natural extension of the pn-QRPA model, developed by Halbleib and Sorensen [Nucl. Phys. A 98, 542 (1967)] to describe the single β decays (A, Z ) → (A, Z ± 1), to the DCE processes. As such, it exhibits several advantages over the pn-QRPA model when used in the evaluation of the double beta decay (DBD) rates. For instance, (i) the extreme sensitivity of the nuclear matrix elements (NMEs) on the model parametrization does not occur; (ii) it allows us to study the NMEs, not only for the ground state in daughter nuclei, as the pn-QRPA model does, but also for all final 0 + and 2 + states, accounting at the same time for their excitation energies and the corresponding DBD Q values; (iii) together with the DBD-NMEs it also provides the energy spectra of Fermi and Gamow-Teller DCE transition strengths, as well as the locations of the corresponding resonances and their sum rules; (iv) the latter are relevant for both the DBD and the DCE reactions, since the underlying nuclear structure is the same; this correlation does not exist within the pn-QRPA model. As an example, detailed numerical calculations are presented for the (A, Z ) → (A, Z + 2) process in 48 Ca → 48 Ti and the (A, Z ) → (A, Z − 2) process in 96 Ru → 96 Mo, involving all final 0 + states and 2 + states.