The advent of machine learning and the availability of big data brought a novel approach for researchers to discover fundamental laws of motion. Computers allow to quickly find underlying physical laws from experimental data, without having in-depth knowledge of the system. Applications are widespread among numerous fields such as physics, chemistry, engineering, biology, climate science and finance. However, the problem is that often experimental data are polluted with noise. This causes symbolic regressors to capture the noise but not the underlying physics. In this work, we leverage the noise reduction properties of non-parametric machine learning, to improve symbolic regression methods. This combination allows for numerically robust differentiation and significantly increases the noise tolerance of common symbolic regressors. A new strategy is exemplified by combining Gaussian processes and the symbolic regression toolbox SINDy for finding governing equations from data. The method balances the quality of the fit and parsimony to avoid overfitting. The method is tested on well-known dynamical systems with varying properties. These will include the Duffing oscillator, as an example of a forced system, the Van der Pol oscillator, as an example of an autonomous system and the Lorenz attractor as an example of a chaotic system.