In this work, the two-center Dirac equation is solved numerically using an extension of an adapted B-spline basis set method previously implemented in relativistic atomic calculations (Fischer, C. F.; Zatsarinny, O. Comput. Phys. Commun. 2009, 180, 879). The robustness of the chosen numerical method, which avoids the appearance of spurious states common in other approaches, allows us to investigate molecular photoionization within a relativistic framework by simply adapting those methods already available in the nonrelativistic case (Brosolo, M.; Decleva, P. Chem. Phys. 1992, 159, 185; Brosolo, M.; Decleva, P.; Lisini, A. Mol. Opt. Phys. 1992, 25, 3345). First, light diatomic molecules (i.e., H₂⁺ and HeH²⁺) are investigated with the purpose of testing the validity and efficiency of the method. Then, a series of one-electron molecular hydrides (i.e., HF⁹⁺, HCl¹⁷⁺ and HI⁵³⁺) is explored by computing the total photoionization cross sections, asymmetry β-parameters and partial phase shifts. The present methodology can be easily extended to treat N-electron molecules following previous approaches in nonrelativistic calculations (Plesiat, E.; Decleva, P.; Martin, F. Phys. Chem. Chem. Phys. 2012, 14, 10853). The inclusion of a second photon can be also accomplished just like in atomic investigations aiming at reproducing pump–probe experiments capable to extract the photoionization time-delays (Vinbladh, J.; Dahlstrom, J. M.; Lindroth, E. Phys. Rev A 2019, 100, 043424; Vinblach, J.; Dahlstrom, J. M.; Lindroth, E. Atoms 2022, 10, 80).