We consider inference for misaligned multivariate functional data that represents the same underlying curve, but where the functional samples have systematic differences in shape. We introduce a class of generally applicable models where warping effects are modelled through non-linear transformation of latent Gaussian variables and systematic shape differences are modelled by Gaussian processes. To model cross-covariance between sample co-ordinates we propose a class of low dimensional cross-covariance structures that are suitable for modelling multivariate functional data. We present a method for doing maximum likelihood estimation in the models and apply the method to three data sets. The first data set is from a motion tracking system where the spatial positions of a large number of body markers are tracked in three dimensions over time. The second data set consists of longitudinal height and weight measurements for Danish boys. The third data set consists of three-dimensional spatial hand paths from a controlled obstacle avoidance experiment. We use the method to estimate the cross-covariance structure and use a classification set-up to demonstrate that the method outperforms state of the art methods for handling misaligned curve data.