We introduce a migration algorithm based on paraxial wave equation that does not use any splitting in the lateral variables. The discretization is first derived in the constant coefficient case by higher order finite differences, then generalized to arbitrarily varying velocities via finite elements. We present a detailed plane wave analysis in a homogeneous medium, and give evidence that numerical dispersion and anisotropy can be controlled. Propagation along depth is done with a higher order method based on a conservative Runge Kutta method. At each step in depth we have to solve a large linear system. This is the most time consuming part of the method. The key to obtaining good performance lies in the use of a Conjugate Gradient like iterative solver. We show the performance of the method with a model example.