We present an adaptation of the fast Fourier factorization method to the simulation of two-dimensional (2D) photonic crystals with a third-order nonlinearity. The algorithm and its performance are detailed and illustrated via the simulation of a Kerr 2D photonic crystal. A change in the transmission spectrum at high intensity is observed. We explain why the change does not reduce to a translation (redshift) but rather consists in a deformation and why one side of the bandgap is more suited to a switching application than the other one.