In this work, Galerkin approximations are developed for a system of first order nonlinear neutral delay differential equations (NDDEs). The NDDEs are converted into an equivalent system of hyperbolic partial differential equations (PDEs) along with the nonlinear boundary constraints. Lagrange multipliers are introduced to enforce the boundary constraints. The explicit expressions for the Lagrange multipliers are derived by exploiting the equivalence of partial derivatives in space and time at a given point on the domain. To illustrate the method, comparisons are made between numerical solution of NDDEs and its Galerkin approximations for different NDDEs. © 2013 American Society of Mechanical Engineers.