A compact numerical method for simulating ultrafast pulse interaction with inhomogeneously broadened multi-level media is reported. We use a low-dispersion pseudospectral scheme with fourth order time stepping for Maxwell's equations, and a weakly coupled operator splitting method for the Bloch equations where inhomogeneous broadening and relaxations are also taken into account. The underlying physics is briefly discussed with emphasis on the formalism used.