While most of the modern seismic imaging methods perform imaging by separating input data into parts (shot gathers), we develop a formulation that is able to incorporate all available data at once while numerically propagating the recorded multidimensional wavefield backward in time. While computationally extensive, this approach has the potential of generating accurate images, free of artifacts associated with conventional approaches. We derive novel high-order partial differential equations in source-receiver-time domain. The fourth order nature of the extrapolation in time has four solutions two of which correspond to the ingoing and outgoing P-waves and reduces to the zero-offset exploding-reflector solutions when the source coincides with the receiver. Using asymptotic approximations, we develop an approach to extrapolating the full prestack wavefield forward or backward in time.