The problem of representing linear systems through combination of positive systems is relevant when signal processing schemes, such as filters, state observers, or control laws, are to be implemented using "positive" technologies, such as Charge Routing Networks and fiber optic filters. This problem, well investigated in the SISO case, can be recasted into the more general problem of Internally Positive Representation (IPR) of systems. This paper presents a methodology for the construction of such IPRs for MIMO systems, based on a suitable convex positive representation of complex vectors and matrices. The stability properties of the IPRs are investigated in depth, achieving the important result that any stable system admits a stable IPR of finite dimension. A main algorithm and three variants, all based on the proposed methodology, are presented for the construction of stable IPRs. All of them are straightforward and are characterized by a very low computational cost. The first and second may require a large state-space dimension to provide a stable IPR, while the third and the fourth are aimed to provide stable IPRs of reduced order.