Novel time domain finite element methods are proposed to numerically solve the system of Maxwell’s equations with a cubic nonlinearity in the spatial 3D case. The effects of linear and nonlinear electric polarization are precisely modeled in this approach. In order to achieve an energy stable discretization at the semi-discrete and the fully discrete levels, a novel technique is developed to handle the discrete nonlinearity, with spatial discretization either using edge and face elements (Nédélec-Raviart-Thomas) or discontinuous spaces and edge elements (Lee-Madsen). In particular, the proposed time discretization scheme is unconditionally stable with respect to the electromagnetic energy and is free of any Courant-Friedrichs-Lewy-type condition. Optimal error estimates are presented at semi-discrete and fully discrete levels for the nonlinear problem. The methods are robust and allow for discretization of complicated geometries and nonlinearities of spatially 3D problems that can be directly derived from the full system of nonlinear Maxwell’s equations.