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 star clusters. Systematic studies dedicated to exploring uncertain parameters in stellar evolution are required to account for the recent observations of the stellar populations. While fitting formulae to stellar tracks in the form of Single Star Evolution (SSE) code remain a popular choice for modelling stellar evolution in population synthesis codes, they are less adaptable to changes in the stellar tracks. Hence, we have developed a Method of Interpolation for Single Star Evolution (METISSE) as an alternative to SSE. It makes use of interpolation between sets of pre-computed stellar tracks to approximate evolution parameters for a population of stars. METISSE is comparable to SSE in performance and can reproduce tracks from different stellar evolution codes quite accurately. In this work, we apply METISSE with detailed stellar tracks computed by the Modules for Experiments in Stellar Astrophysics (MESA), Bonn Evolutionary Code (BEC), to study the impact of uncertainties in stellar evolution on a population 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, including remnant masses and maximal radial expansion. The differences in the predictions of different stellar models can help us account for the present day observations of stellar populations.