Non-Born-Oppenheimer variational calculations employing explicitly correlated Gaussian basis functions have been performed for the ground states of the beryllium monohydride molecule (BeH) and its ion (BeH+), as well as for the beryllium atom (Be) and its ion (Be+). An approach based on the analytical energy gradient calculated with respect to the Gaussian exponential parameters was employed. The calculated energies were used to determine the ionization potential of BeH and the dissociation energies of BeH and BeH+. Also, the generated wave functions were used to compute various expectation values, such as the average interparticle distances and the nucleus-nucleus correlation functions.