We revisit the non-relativistic limit of the Dirac equation in finite scalar and vector potentials and derive a Schrödinger-like equation that retains leading spin–relativistic corrections in closed form. For general central potentials, we cast the radial equation into a quadratic eigenvalue problem (QEP) using a finite-difference discretization method and develop an open-source solver to address it. We study Coulomb, harmonic oscillator, Woods–Saxon, and Yukawa potentials. We further obtain first-order energy and wavefunction corrections for the three-dimensional isotropic harmonic oscillator and Coulomb potentials via perturbation theory. This framework provides a practical bridge between non-relativistic and fully relativistic treatments, enabling accurate quantification of relativistic effects without the computational cost of full four-component calculations.