SummaryThis paper proposes a new numerical approach for finding the solution of linear time-delay control systems with a quadratic performance index using new hybrid functions. This method is based on a hybrid of block-pulse functions and biorthogonal multiwavelets that consist of cubic Hermite splines on the primal side. The excellent properties of the hybrid functions together with the operational matrices of integration, product, and delay are presented. Using these matrices, the solution of the optimal control of delay systems is reduced to the solution of algebraic equations. Because of the sparsity nature of these matrices, this method is computationally very attractive and reduces CPU time and computer memory. In order to save the memory requirement and computation time, a threshold procedure is applied to obtain algebraic equations. The effectiveness of the method and the accuracy of the solution are shown in comparison with some other methods by illustrative examples.