  The mass spectrum of heavy pseudoscalar mesons, described as quark-antiquark
bound systems, is considered within the Bethe-Salpeter formalism with momentum
dependent masses of the constituents. This dependence is found by solving the
Schwinger-Dyson equation for quark propagators in rainbow-ladder approximation.
Such an approximation is known to provide both a fast convergence of numerical
methods and accurate results for lightest mesons. However, as the meson mass
increases, the method becomes less stable and special attention must be devoted
to details of numerical means of solving the corresponding equations. We focus
on the pseudoscalar sector and show that our numerical scheme describes fairly
accurately the $\pi$, $K$, $D$, $D_s$ and $\eta_c$ ground states. Excited
states are considered as well. Our calculations are directly related to the
future physics programme at FAIR.
