Commonly, variance-based global sensitivity analysis methods are popular and applicable to quantify the impact of a set of input variables on output response. However, for many engineering practical problems, the output response is not single but multiple, which makes some traditional sensitivity analysis methods difficult or unsuitable. Therefore, a novel global sensitivity analysis method is presented to evaluate the importance of multi-input variables to multi-output responses. First, assume that a multi-input multi-output system (MIMOS) includes [Formula: see text] variables and [Formula: see text] responses. A set of summatory functions [Formula: see text] and [Formula: see text] are constructed by the addition and subtraction of any two response functions. Naturally, each response function is represented using a set of summatory function. Subsequently, the summatory functions [Formula: see text] and [Formula: see text] are further decomposed based on the high dimensional model representation (HDMR), respectively. Due to the orthogonality of all the decomposed function sub-terms, the variance and covariance of each response function can be represented using the partial variances of all the decomposed function sub-terms on the corresponding summatory functions, respectively. The total fluctuation of MIMOS is calculated by the sum of the variances and covariances on all the response functions. Further, the fluctuation is represented as the sum of the total partial variances for all the [Formula: see text]-order function sub-terms, and the total partial variance is the sum of [Formula: see text] partial variances for the corresponding [Formula: see text]-order function sub-terms. Then, the function sensitivity index (FSI) [Formula: see text] for s-order function sub-terms is defined by the ratio of the total partial variance and total fluctuation, which includes first-order, second-order, and high-order FSI. The variable sensitivity index [Formula: see text] of variable [Formula: see text] is calculated by the sum of all the FSIs including the contribution of variable [Formula: see text]. Finally, numerical example and engineering application are employed to demonstrate the accuracy and practicality of the presented global sensitivity analysis method for MIMOS.