The numerical approximation of the solution of the time-dependent Schrödinger equation arising in ultrafast laser dynamics is discussed. The linear electronic Schrödinger equation is reduced to a computationally tractable, lower dimensional system of nonlinear partial differential equations by the \emph{multi-configuration time-dependent Hartree-Fock method (MCTDHF)}. This method uses an ansatz for the wave function on a nonlinear manifold, taking into account the antisymmetry inherent in the model to reduce the dimension of the solution space. We show that the ansatz used is consistent with the original equation and the solution can be represented exactly in principle using the solutions of the nonlinear PDEs associated with MCTDHF. For the practical solution of the MCTDHF equations, several numerical techniques are discussed, and it is demonstrated that physically relevant problems can be solved satisfactorily.