A new method for identification of systems of arbitrary real order based on numerical solution of systems of nonlinear fractional order differential equations (FODEs) and orthogonal dis-tance fitting is presented. The main idea is to fit experimental or measured data using a solution of a system of fractional differential equations. The parameters of these equations, including the orders of derivatives, are subject to optimization process, where the criterion of optimization is the minimal sum of orthogonal distances of the data points from the fitting line. Once the minimal sum is found, the identified parameters are considered as optimal. The so called orthogonal distance fitting, known also under the names of total least squares or orthogonal regression is naturally used in the fitting criterion, since it is the most suitable tool for fitting lines and surfaces in multidimensional space. The examples illustrating the methods are presented in 2-dimensional and 3-dimensional problems.