  We have developed a formalism, based on the Fourier-Bessel expansion, that
facilitates the evaluation of matrix elements involving nucleon recoil
operators, such as appear in serveral exotic forms of neutrinoless double beta
decay ($\beta\beta_{0\nu}$). The method is illustrated by applying it to the
``charged'' majoron model, which is one of the few that can hope to produce an
observable effect. From our numerical computations within the QRPA performed
for $^{76}Ge$, $^{82}Se$, $^{100} Mo$, $^{128}Te$ and $^{150}Nd$ nuclei, we
test the validity of approximations made in earlier work to simplify the new
matrix elements, showing that they are accurate to within 15%. Our new method
is also suitable for computing other previously unevaluated $\beta\beta_{0\nu}$
nuclear matrix elements.
