In the era of advanced electromagnetic and gravitational wave detectors, it has become increasingly important to effectively combine and study the impact of stellar evolution on binaries and dynamical systems of stars. Systematic studies dedicated to exploring uncertain parameters in stellar evolution are required to account for the recent observations of the stellar populations. We present a new approach to the commonly used single-star evolution (sse) fitting formulae, one that is more adaptable: method of interpolation for single star evolution (metisse). It makes use of interpolation between sets of pre-computed stellar tracks to approximate evolution parameters for a population of stars. We have used metisse with detailed stellar tracks computed by the modules for experiments in stellar astrophysics (mesa), the bonn evolutionary code (bec), and the Cambridge stars code. metisse better reproduces stellar tracks computed using the stars code compared to sse, and is on average three times faster. Using stellar tracks computed with mesa and bec, we apply metisse to explore the differences in the remnant masses, the maximum radial expansion, and the main-sequence lifetime of massive stars. We find that different physical ingredients used in the evolution of stars, such as the treatment of radiation-dominated envelopes, can impact their evolutionary outcome. For stars in the mass range 9–100 M_⊙, the predictions of remnant masses can vary by up to 20 M_⊙, while the maximum radial expansion achieved by a star can differ by an order of magnitude between different stellar models.