In this thesis, an efficient numerical method for solving nonlinear optimal control of multi-delay systems with inequality constraints on the state and control variables is presented. Many phenomena arising in various fields of science and engineering can be described by delay differential equations. Typical examples are chemical processes, transmission lines, biological models, population growth, economics and finance and communication networks. So far several computational techniques have been devoted to the numerical solution of delay differential equations. Optimal control of multi-delay systems is one of the most challenging problems i control theory. It is known that except for simple cases, it is either difficult or impossible to analytically solve a nonlinear optimal control problem involving constraints on the state and control variables. As a result, an effective numerical method has to be developed. Numerical methods for solving nonlinear optimal control problems are typically described into two categories: direct methods and indirect methods. Historically, many early numerical methods were based on finding solutions to satisfy a set of necessary optimally conditions resulting from Pontryagin’s maximum principle. These methods are collectively called indirect methods. There are many successful implementation of indirect methods in the literature. Although indirect methods enjoy some nice properties, they suffer from many drawbacks. For instance, the boundary value problem resulting from the necessary optimally conditions is extremely sensitive to initial guesses. In addition, these necessary conditions must be explicitly derived. Over the past two decades, an alternative approach based on discrete approximations has gained wide popularity. The essential idea of this method is to discretize the optimal control problem and solve the resulting finite-dimensional optimization problem. The simplicity of direct methods belies a wide range of deeply theoretical issues that lie at the intersection of approximation theory, control theory and optimization. The method implemented in this thesis is based on direct approach using a hybrid of block-pulse function and Legendre polynomials using the well-known Legendre-Gauss. A penalty function method is employed to convert the inequality constraints on the state and control variables into equality constraints. A simple strategy is devised for obtaining the required number of subintervals. The suitable value of N, the order of block-pulse functions enables one to accurately adjust the exact locations of the switching points that occur in the control input. After choosing the value of N, we can improve the accuracy of the solution by increasing the value of M, the degree of the Legendre polynomials. It should be noted that the analytical solution associated with a delayed optimal control problem is a piecewise smooth function. For this reason, smooth basis functions such as Legendre polynomials are not able to provide a satisfactory approximation for time-delay systems. The nice properties of the hybrid functions are used to transform the delayed optimal control problem into a parameter optimization problem whose solution is much easier than the original one. Several numerical examples are investigated to show the efficiency and accuracy of the proposed numerical method. The superiority of the hybrid functions is demonstrated over the Legendre polynomials.