Multivariate functional data from a complex system are naturally high-dimensional and have 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 paper presents a dynamic subspace learning method for multivariate functional data modeling. In particular, we consider 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 (FISTA), and the features of every subspace can be extracted using the smooth multi-channel functional PCA. Numerical studies together with case studies demonstrate the efficiency and applicability of the proposed methodology.