A numerical method to determine the stability of delay differential equations (DDEs) with time periodic coefficients is proposed. The DDE is converted into an equivalent partial differential equation (PDE) with a time periodic boundary condition (BC). The PDE, along with its BC, is then converted into a system of ordinary differential equations (ODEs) with time periodic coefficients using the Galerkin least squares approach. In the Galerkin approach, shifted Legendre polynomials are used as basis functions, allowing us to obtain explicit expressions for the approximate system of ODEs. We analyze the stability of the discretized ODEs, which represent an approximate model of the DDEs, using Floquet theory. We use numerical examples to show that the stability charts obtained with our method are in excellent agreement with those existing in the literature and those obtained from direct numerical simulation.