In this paper, we develop a method for accurately modeling ultrafast molecular dynamics and pump-probe spectroscopy using classical trajectory simulations. The approach is based on a semiclassical limit of the Liouville formulation of quantum mechanics. Expressions for the nonstationary classical phase space probability density created by an ultrashort laser pulse on an excited electronic state, and the observable fluorescence signal resulting from a pump-probe experiment, are derived in the weak-held limit using perturbation theory. By introducing additional approximations, these expressions are cast in a form that can be directly implemented using classical trajectory integration and ensemble averaging. The method is tested against multisurface time-dependent quantum mechanical wave packet calculations for a one-dimensional model system representing I-2 photodissociation-recombination in a static Ar lattice. Nearly quantitative agreement between the exact calculations and the trajectory-based method is obtained. Although demonstrated for a one-dimensional system, the method is easily incorporated in conventional molecular dynamics programs, allowing efficient treatment of many-body ultrafast dynamics and spectroscopy.