  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
solutions, parallelization, and, in some cases, a lower 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 $\Lambda \mathrm{CDM}$, the
Chevallier-Polarski-Linder parametric dark energy model, a quintessence model
with an exponential potential, and the Hu-Sawicki $f(R)$ model. We use 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, type Ia supernovae
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 estimate the
error of the solutions that the trained networks provide by comparing them with
the analytical solution when there is one, or to a high-precision numerical
solution when there is not. Through these estimations we find that the error of
the solutions is at most $\sim1\%$ in the region of the parameter space that
concerns the $95\%$ confidence regions that we find using the data, for all
models and all statistical analyses performed in this work.
