In the last few years, there has been significant progress in the development
of machine learning methods tailored to astrophysics and cosmology. Among the
various methods that have been developed, there is one that allows to obtain a
bundle of solutions of differential systems without the need of using
traditional numerical solvers. We have recently applied this to the
cosmological scenario and showed that in some cases the computational times of
the inference process can be reduced. In this paper, we present an improvement
to the neural network bundle method that results in a significant reduction of
the computational times of the statistical analysis. The novelty of the method
consists in the use of the neural network bundle method to calculate the
luminosity distance of type Ia supernovae, which is usually computed through an
integral with numerical methods. In this work, we have applied this improvement
to the Starobinsky f(R) model, which is more difficult to integrate than the
f(R) models analyzed in our previous work. We performed a statistical
analysis with data from type Ia supernovae of the Pantheon+ compilation and
cosmic chronometers to estimate the values of the free parameters of the
Starobinsky model. We show that the statistical analyses carried out with our
new method require lower computational times than the ones performed with both
the numerical and the neural network method from our previous work. This
reduction in time is more significant in the case of a difficult computational
problem such as the one we address in this work.Comment: 11 pages, 3 figures, 2 tables, to be submitted to PR