Multivariate functional data from a complex system are naturally high-dimensional and have a complex cross-correlation structure. The complexity of data structure can be observed as that (1) some functions are strongly correlated with similar features, while some others may have almost no cross-correlations with quite diverse features; and (2) the cross-correlation structure may also change over time due to the system evolution. With this regard, this article presents a dynamic subspace learning method for multivariate functional data modeling. In particular, we consider that different functions come from different subspaces, and only functions of the same subspace have cross-correlations with each other. The subspaces can be automatically formulated and learned by reformatting the problem as a sparse regression. By allowing but regularizing the regression change over time, we can describe the cross-correlation dynamics. The model can be efficiently estimated by the fast iterative shrinkage-thresholding algorithm, and the features of each subspace can be extracted using the smooth multi-channel functional principal component analysis. Some theoretical properties of the model are presented. Numerical studies, together with case studies, demonstrate the efficiency and applicability of the proposed methodology.