Understanding treatment heterogeneity plays a key role in many contemporary applications arising from different areas. Although there is a growing literature on subgroup analysis based on the heterogeneous univariate regression model, little work has been done for the heterogeneous multi-response regression model. In contrast to the existing methods which are based on subject-specific intercepts, we introduce a heterogeneous multi-response regression method which allows the coefficients for treatment variables to be subject-dependent with unknown grouping information, thus applicable to a wider range of situations. Moreover, we provide an efficient algorithm based on concave pairwise fusion penalization and establish the oracle property of the proposed estimator. The effectiveness of the suggested method is demonstrated through simulation examples and an empirical study.