This thesis is devoted to the numerical computation of the time-independent Dirac equation (TIDE) and the time-dependent Dirac equation (TDDE) in the prolate spheroidal coordinates. Analytical and numerical techniques including Galerkin methods, Min-max principle, and Rayleigh-Ritz methods combined with atomically balanced basis are presented to solve the Dirac equation without spectral pollution. These numerical methods are used to compute the discrete spectrum of the Dirac operator in two-center Coulomb problems for molecules $\mbox{H}_2^+$ and $\mbox{Th}_2^{179+}$. High order B-spline basis functions are used to obtain accurate results. As excepted, the numerical results do not show any spurious state.