Correlation functions of Yang-Mills theory in the Landau gauge are calculated from their equations of motion. The employed setup is completely parameter-free and leads, within errors, to good quantitative agreement with corresponding lattice results for the ghost and gluon propagators as well as the ghost-gluon and three-gluon vertices. Also the four-gluon vertex is calculated. The present setup allows for the first time for a unique subtraction of quadratic divergences in the gluon propagator Dyson-Schwinger equation. Thus, there are no ambiguities which can arise due to the use of models or auxiliary workarounds. In addition, several self-tests of the results are described that allow assessing the truncation error in a self-consistent way. This enables a new perspective on how to identify limitations of the present setup and develop future improvements.