An algebraic solution of the Dirac equation is reinvestigated. Slater-type spinor orbitals and their corresponding system of differential equations are defined in two- and four-component formalism. They describe the radial function in components of the wave function of the Dirac equation solution to high accuracy. They constitute the matrix elements arising in a generalized eigenvalue equation. These terms are evaluated through prolate spheroidal coordinates. The corresponding integrals are calculated by the numerical global-adaptive method taking into account the Gauss-Kronrod numerical integration extension. Sample calculations are performed using flexible basis sets generated with both signs of the relativistic angular momentum quantum number κ. Applications to one-electron atoms and diatomics are detailed. Variationally optimum values for orbital parameters are obtained at given nuclear separation. Methods discussed in this work are capable of yielding highly accurate relativistic two-center integrals for all ranges of orbital parameters. This work provides an efficient way to overcome the problems that arise in relativistic calculations.