The field of machine learning has drawn increasing interest from various other fields due to the success of its methods at solving a plethora of different problems. An application of these has been to train artificial neural networks to solve differential equations without the need of a numerical solver. This particular application offers an alternative to conventional numerical methods, with advantages such as lower memory required to store the solutions, parallelization, and in some cases less overall computational cost than its numerical counterparts. In this work, we train artificial neural networks to represent a bundle of solutions of the differential equations that govern the background dynamics of the Universe for four different models. The models we have chosen are ΛCDM, the Chevallier-Polarski-Linder parametric dark energy model, a quintessence model with an exponential potential, and the Hu-Sawicki f (R) model. We used the solutions that the networks provide to perform statistical analyses to estimate the values of each model's parameters with observational data; namely, estimates of the Hubble parameter from Cosmic Chronometers, the Supernovae type Ia data from the Pantheon compilation, and measurements from Baryon Acoustic Oscillations. The results we obtain for all models match similar estimations done in the literature using numerical solvers. In addition, we estimated the error of the solutions by comparing them to the analytical solution when there is one, or to a high-precision numerical solution when there is not. Through those estimations we found that the error of the solutions was at most ∼ 1% in the region of the parameter space that concerns the 95% confidence regions that we found using the data, for all models and all statistical analyses performed in this work.