We present an algorithm for data-driven reachability analysis that estimates finite-horizon forward reachable sets for general nonlinear systems using level sets of a certain class of polynomials known as Christoffel functions. The level sets of Christoffel functions are known empirically to provide good approximations to the support of probability distributions: the algorithm uses this property for reachability analysis by solving a probabilistic relaxation of the reachable set computation problem. We also provide a guarantee that the output of the algorithm is an accurate reachable set approximation in a probabilistic sense, provided that a certain sample size is attained. We also investigate three numerical examples to demonstrate the algorithms capabilities, such as providing non-convex reachable set approximations and detecting holes in the reachable set.