The numerical approximation of the solution of the time-dependent Schrödinger equation arising in ultrafast laser dynamics is discussed. The linear Schrödinger equation is reduced to a computationally tractable, lower dimensional system of nonlinear partial differential equations by the multi-configuration time-dependent Hartree-Fock method. This method serves to approximate the original wave function on a nonlinear manifold, using the antisymmetry inherent in the model to significantly reduce the dimension of the solution space. For the solution of the resulting systems of PDEs, several numerical techniques are compared. Space discretization using the pseudospectral method turns out to be superior to finite difference approximations. For time integration, the range of applicability and computational efficiency of high-order Runge-Kutta methods are compared with variational splitting, a method recently proposed for quantum molecular dynamics.