Using the operational matrix of an orthogonal function to perform integration for solving, identifying and optimizing a linear dynamic system has several advantages: (1) the method is computer oriented, thus solving higher order differential equation becomes a matter of dimension increasing; (2) the solution is a multiresolution type; (3) the answer is convergent, even the size of increment is very large. The traditional methods of deriving the operational matrix is much involved and not unified, this paper presents a new unified approach to deriving the operational matrices of orthogonal functions. We apply it first to the derivation of the operational matrices of the square wave group which consist of (i) the block pulse function, (ii) the Walsh function and (iii) the Haar wavelet function, then to the sinusoidal group which includes (i) the discrete Fourier transform, (ii) the discrete cosine transform and (iii) the discrete Hartley transform. Finally, we use the operational matric...