Machine learning strategies can be efficiently used to accelerate the exploration of the design space or flight envelope of highly flexible aeroelastic systems. In this paper, we explore the use of interpolation between parametric state-space realizations to, with few true systems sampled in the parameter space, produce with adequate accuracy a state-space model anywhere in the parameter space. The location of the sampling points is shown to be decisive thus the selection of these points takes the focus in this work. Several approaches are explored, putting emphasis on adaptive schemes that locate the optimal points in the parameter space that are needed to capture the changing system dynamics. Since the evaluation of the true system is costly, optimization techniques based on statistical surrogate models are sought, which need to be trained but are effective in locating the best locations to use as sampling data. A novel method inspired by Bayesian optimization is used to make the most out of a limited number of known state-spaces by taking different combinations as training and testing data of the statistical surrogate, leading to not only an accurate interpolation framework but also to a reduction of 50% of the number of costly full system evaluations compared to a standard Bayesian optimization set-up. These methods are demonstrated on the Pazy wing, a very flexible wing with a complex stability envelope, whereby we produce a very accurate representation of the flutter envelope at a reduced computational cost.