One of the benefits of using Gaussian Process models is the availability of uncertainty quantification. However, when learning continuous dynamical systems obtaining trajectories requires repeatedly mapping uncertain inputs through the learned nonlinear function, which is generally non-tractable. As sampling-based approaches are computationally expensive, we consider approximations of the output and trajectory distribution. We show that existing approaches make an incorrect implicit independence assumption and underestimate the model-induced uncertainty. We propose a piecewise linear approximation of the GP model and a numerical solver for efficient uncertainty estimates matching sampling-based methods at a lower computational cost.