Potential energy surfaces of the hydrogen molecular ion H$_2^+$ in the Born-Oppenheimer approximation are computed by means of the Riccati-Pade method (RPM). The convergence properties of the method are analyzed for different states. The equilibrium internuclear distance, as well as the corresponding electronic plus nuclear energy, and the associated separation constants, are computed to 40 digits of accuracy for several bound states. For the ground state the same parameters are computed with more than 100 digits of accuracy. Additional benchmark values of the electronic energy at different internuclear distances are given for several additional states. The software implementation of the RPM is given under a free software license. The results obtained in the present work are the most accurate available so far, and further additional benchmarks are made available through the software provided.