Bifurcation theory tells us that dynamics can be sensitive to the choice of model. This suggests a potential inherent instability in the following commonly used data analysis pipeline: use time series data to identify/fit an explicit model, e.g., a differential equation or map, and then use the model to predict the dynamics. One expects that this potential instability is magnified by sparse data or noisy data. We propose a novel method, combining Conley theory and Gaussian Process surrogate modeling with uncertainty quantification, through which it is possible to characterize local and global dynamics, e.g., existence of fixed points, periodic orbits, connecting orbits, bistability, and chaotic dynamics, and to provide lower bounds on the confidence that this characterization of the dynamics is correct. Furthermore, numerical experiments indicate that it is possible to identify nontrivial dynamics with high confidence with surprisingly small data sets.