Synchrophasor data provide unprecedented opportunities for inferring power system dynamics, such as estimating voltage angles, frequencies, and accelerations along with power injection at all buses. Aligned to this goal, this work puts forth a novel framework for learning dynamics after small-signal disturbances by leveraging Gaussian processes (GPs). We extend results on learning of a linear time-invariant system using GPs to the multi-input multi-output setup. This is accomplished by decomposing power system dynamics into a set of single-input single-output linear systems with narrow frequency pass bands. The proposed learning technique captures time derivatives in continuous time, accommodates data streams sampled at different rates, and can cope with missing data and heterogeneous levels of accuracy. While Kalman filter-based approaches require knowing all system inputs, the proposed framework handles readings of system inputs, outputs, their derivatives, and combinations thereof collected from an arbitrary subset of buses. Relying on minimal system information, it further provides uncertainty quantification in addition to point estimates of system dynamics. Numerical tests verify that this technique can infer dynamics at non-metered buses, impute and predict synchrophasors, and locate faults under linear and non-linear system models under ambient and fault disturbances.