The problems of systems identification, analysis and optimal control have been recently studied using orthogonal functions. The specific orthogonal functions used up to now are the Walsh, the block-pulse, the Laguerre, the Legendre, Haar and many other functions. In the present paper, several operational matrices for integration and differentiation are studied. we introduce the Kronecker convolution product and expanded to the Riemann–Liouville fractional integral of matrices. For some applications, it is often not necessary to compute exact solutions, approximate solutions are sufficient because sometimes computational efforts rapidly increase with the size of matrix functions. Our method is extended to find the exact and approximate solutions of the general system matrix convolution differential equations, the way exists which transform the coupled matrix differential equations into forms for which solutions may be readily computed. Finally, several systems are solved by the new and other approaches and illustrative examples are also considered.