This paper presents a novel approach for propagating uncertainties in dynamical systems building on high-order Taylor expansions of the flow and moment-generating functions (MGFs). Unlike prior methods that focus on Gaussian distributions, our approach leverages the relationship between MGFs and distribution moments to extend high-order uncertainty propagation techniques to non-Gaussian scenarios. This significantly broadens the applicability of these methods to a wider range of problems and uncertainty types. High-order moment computations are performed one-off and symbolically, reducing the computational burden of the technique to the calculation of Taylor series coefficients around a nominal trajectory, achieved by efficiently integrating the system's variational equations. Furthermore, the use of the proposed approach in combination with event transition tensors, allows for accurate propagation of uncertainties at specific events, such as the landing surface of a celestial body, the crossing of a predefined Poincaré section, or the trigger of an arbitrary event during the propagation. Via numerical simulations we demonstrate the effectiveness of our method in various astrodynamics applications, including the unperturbed and perturbed two-body problem, and the circular restricted three-body problem, showing that it accurately propagates non-Gaussian uncertainties both at future times and at event manifolds.