Time-delay systems have received much attention in the past few decades because time delays are frequently encountered in many practical systems and various fields of engineering and science, such as aerospace engineering, robotics, physics, communication networks, chemical processes, traortation systems, transmission lines, bilogical models, population growth, economics and finance, climate models and power systems. The presence of delay makes analysis and control design much more complicated. Therefore, time-delay systems are very important to many investigators for their control, stability, and optimization. The application of Pontryagin’s maximum principle to the optimization of control systems with time delays as outlined by Kharatishvili results in a system of coupled two-point boundary value problem involving both delayed and advanced terms whose exact solution, except in highly special cases, is very difficult. Therefore, the main object of all computational aspects of optimal time-delay systems has been to devise a methodology to avoid the solution of the mentioned two-point boundary value problem. Orthogonal functions have been extensively used to solve various problems of dynamic systems. The essential idea of this technique is that it reduces these problems to those of solving a system of algebraic equations, thus greatly simplifying the problem. Signals frequently have mixed features of continuity and jumps. These signals are continuous over certain segments of time with discontinuities or jump occurring at the transitions of the segments. In such situations, neither the continuous basis functions nor piecewise constant basis functions taken alone would form an efficient basis in the representation of such signals. In general, the computed response of the delay systems via continuous orthogonal functions or piecewise constant basis functions is not in good agreement with the exact response of the system. In this thesis, a composite Chebyshev finite difference method is applied for finding the solution of optimal control of time-delay systems with a quadratic performance index. The proposed approximation method is based on a hybrid of block-pulse functions and Chebyshev polynomials using the Chebyshev-Gauss-Lobatto points. Because the set of block-pulse functions and Chebyshev polynomials are both complete and orthogonal, the set of hybrid functions is a complete orthogonal set. An essential feature of the proposed hybrid functions is the good representation of smooth and especially piecewise smooth functions. It should be mentioned that the exact solution of delayed optimal control problems cannot be obtained by smooth basis functions. The suggested method is an extension of the Chebyshev finite difference scheme. The proposed method can be regarded as a nonuniform finite difference scheme. The convergence of presented procedure is discussed. The nice properties of hybrid functio are used to convert the optimal control problem into a mathematical programming problem whose solution is much easier than the original one. A wide variety of linear optimal control problems with time delay are included to demonstrate the validity and the applicability of the technique. The method is easy to implement and provides very accurate results. The approximation method developed in this thesis can be applied to linear optimal control problems with multiple time-delays. However, some modifications are required.