We present a method for ?nding the solution of linear multi-delay systems (MDS) by using the hybrid of block-pulse functions and Bernoulli polynomials. In this approach, operational matrices of integration, delay and product are utilized to reduce the solution of multi-delay systems to the solution of a linear system of algebraic equations whose solution is much more easier than the original one. The available sets of orthogonal functions can be divided into three We use the hybrid functions consisting of the combination of the block-pulse functions and the Bernoulli polynomials to solve linear multi-delay systems. The presented method is based upon reducing the multi-delay systems into a set of algebraic equations. We ?rst describe the basic properties of the proposed hybrid functions required for our subsequent development. The excellent properties of the hybrid functions together with the operational matrices of integration and delay are used to convert the multi-delay system into a system of linear algebraic equations. The mentioned operational matrices are sparce. The sparsity of these matrices makes the method computationally attractive at the same time keeping the accuracy of the solution. Illustrative examples are included to demonstrate the validity and applicability of the presented approximation scheme. The simulation results indicate that the method developed in this thesis is e?ective.