Recent improvements in helioseismological data have necessitated an increased accuracy of solar models. The precision of the numerical procedure used to construct a solar model should exceed the highest observational accuracy. Only when this condition is met may differences between models and observations be ascribed to physical assumptions. In this paper 1-M_solar models in various main-sequence phases and solar models that have been calculated by two completely different numerical methods, applying identical physics, are compared. We show that, to achieve high-precision models, a multiple shooting method is better suited than a differencing scheme. The application of such integrator methods to solar model codes should therefore allow a deeper insight into the physical interpretation of helioseismological data.