  We present a new method to solve the nuclear density functional theory (DFT)
equations using a two-center harmonic oscillator basis, incorporating pairing
and Coulomb interactions. The goal is to efficiently determine the fission and
fusion configurations in nuclei. The Coulomb exchange terms are evaluated
exactly, allowing for a novel approach to neck formation without the Slater
approximation, which has been commonly used in space coordinate-based
approaches. The new method has been implemented in the code {\sc hfodd},
enabling direct comparison with standard one-center solutions. This first paper
focuses on deriving and implementing a methodology based on stable, precise,
and exact applications of harmonic oscillator bases for the two fragments,
which can either overlap or be separated by arbitrarily large distances. The
implementation is tested on two proof-of-principle examples using light nuclei,
specifically $^8$Be and $^{24}$Mg.
