Uncertainty propagation in non-linear dynamical systems has become a key problem in various fields including control theory and machine learning. In this work, we focus on discrete-time non-linear stochastic dynamical systems. We present a novel approach to approximate the distribution of the system over a given finite time horizon with a mixture of distributions. The key novelty of our approach is that it not only provides tractable approximations for the distribution of a nonlinear stochastic system but also comes with formal guarantees of correctness. In particular, we consider the Total Variation (TV) distance to quantify the distance between two distributions and derive an upper bound on the TV between the distribution of the original system and the approximating mixture distribution derived from our framework. We show that in various cases of interest, including in the case of Gaussian noise, the resulting bound can be efficiently computed in closed form. This allows us to quantify the correctness of the approximation and to optimize the parameters of the resulting mixture distribution to minimize such distance. The effectiveness of our approach is illustrated on several benchmarks from the control community.@en