Many conventional physical and engineering phenomena have been identified to be well expressed by making use of the fractional order partial differential equations (FOPDEs). For that reason, a proficient and stable numerical method is needed to find the approximate solution of FOPDEs. This article is designed to develop the numerical scheme able to find the approximate solution of generalized fractional order coupled systems (FOCSs) with mixed partial derivative terms of fractional order. Our main objective in this article is the development of a new operational matrix for fractional mixed partial derivatives based on the orthogonal shifted Legendre polynomials (SLPs). The fractional derivatives are considered herein in the sense of Caputo. The proposed method has the advantage to reduce the considered problems to a system of algebraic equations which are simple in handling by any computational software. Being easily solvable, the associated algebraic system leads to finding the solution of the problem. Some examples are included to demonstrate the accuracy and validity of the proposed method.