Coupled Monte Carlo burnup codes aim to evaluate the time evolution of different parameters, such as nuclide densities, for accurate modeling of the different reactor designs and associated fuel cycles. Recently a major deficiency in numerical stability of existing Monte Carlo coupling schemes was identified. Alternative, stable coupling schemes were derived, implemented and verified. These methods are iterative and rely on either the end- or middle-of-step (MOS) reaction rates to evaluate the end-of-step (EOS) nuclide densities. Here, we demonstrate that applying the EOS methods for realistic problems may lead to highly inaccurate results. Considerable improvement can be made by adopting MOS method but the accuracy may still be insufficient. The solution proposed in this work relies on the substep method that allows reducing the time discretization errors. The proposed and tested substep method also assumes that the reaction rates are linear functions of the logarithm of the nuclide densities. The method was implemented in BGCore code and subsequently used to perform a series of test case calculations. The results demonstrate that better accuracy and hence efficiency can be achieved with negligible additional computational burden.