  With the first direct detections of gravitational waves (GWs) from the
coalescence of compact binaries observed by the advanced LIGO and VIRGO
interferometers, the era of GW astronomy has begun. Whilst there is strong
evidence that the observed GWs are connected to the merger of two black holes
(BH) or two neutron stars (NS), future detections may present a less consistent
picture. Indeed, the possibility that the observed GW signal was created by a
merger of exotic compact objects (ECOs) such as boson stars (BS) or axion stars
(AS) has not yet been fully excluded. For a detailed understanding of the late
stages of the coalescence full 3D numerical relativity simulations are
essential. In this paper, we extend the infrastructure of the numerical
relativity code BAM, to permit the simultaneous simulation of baryonic matter
with bosonic scalar fields, thus enabling the study of BS-BS, BS-NS, and BS-BH
mergers. We present a large number of single star evolutions to test the newly
implemented routines, and to quantify the numerical challenges of such
simulations, which we find to partially differ from the default NS case. We
also compare head-on BS-BS simulations with independent numerical relativity
codes, namely the SpEC and the GRChombo codes, and find good general agreement.
Finally, we present what are, to the best of our knowledge, the first full NR
simulations of BS-NS mergers, a first step towards identifying the hallmarks of
BS-NS interactions in the strong gravity regime, as well as possible GW and
electromagnetic observables.
