The study of r-process elements, such as barium (Ba) and strontium (Sr), is key to understanding the formation of heavy elements in the Universe. More precisely, the hyperfine structure (HFS) of their atomic levels is commonly used by astrophysicists to determine abundances via spectral line modeling. Using the Multiconfiguration Dirac–Hartree–Fock (MCDHF) method as implemented in the General Relativistic Atomic Structure Package (GRASP) code, the magnetic dipole and electric quadrupole hyperfine structure constants were determined for the first excited states of Ba II isotopes, as well as for Ba I and Sr II isotopes to monitor the robustness of the developed model. New code developments, such as the use of natural orbitals, the addition of polarization effects and the use of Configuration State Function Generators (CSFGs) as implemented in GRASPG, were tested for these heavy elements. The developed strategy allowed us to achieve encouraging results, with little disagreement with experimental data for all studied level except ²D_5/2 in the first Ba II isotope. This disagreement was observed in another Ba II isotope as well as in Sr II. However, it emerged that the adopted strategy could not describe all the physics crucial for the Ba I states studied, with disagreements reaching up to 70%. This limitation necessitates the introduction of more intensive work employing a multireference (MR) approach to describe configuration mixing. Such efforts are expected to yield improved agreement with experimental values and outcomes from other theoretical computations. Nevertheless, the adopted strategy continues to serve as a valuable benchmark for comparing the MCDHF method with other theoretical approaches.