According to the method of calculating maximal Lyapunov exponent (MLE) from univariate small data sets, an extended method based on multivariate time series is proposed. The extended method can search out optimal reconstructing parameters to meet the requirement of the original method for reconstructing multivariate phase space, and the method can compute the MLE by making use of the optimal reconstructed multivariate phase space. The method is tested by coupled non-identical chaotic Rssler, coupled chaotic Rssler and hyper chaotic Rssler. The test results show that the extend method is efficient, and the computing results of MLE based on multivariate are much closer to the theoretical values than the results of univariate even when the data sets of each time series become small.