Bayesian model selection using Bayes factors requires the computation of marginal likelihoods. Nowadays, the marginal likelihoods are often computed using thermodynamic integration for power posteriors, which relies on numerical integration methods. The commonly used integration methods however neither control the integration accuracy nor exploit the available function evaluations efficiently. In this manuscript we introduce an adaptive method for calculating marginal likelihoods which relies on Simpson’s rule. The proposed method is evaluated on an analytically tractable academic example as well as two high-dimensional models possessing up to 800 parameters. The high-dimensional models describe the protein degradation in a large population of fibroblast cells. Our analysis reveals that the proposed adaptive method shows improved performance over existing approaches for simple problems and furthermore allows for the efficient study of high-dimensional problems.