We proposed a technique for the fast solution of ray tracing in three-dimensional laterally heterogeneous media. Analytical solutions for both ray tracing and paraxial ray tracing for a medium with an arbitrary gradient of the square of slowness are used to develop a systematic procedure to construct exact solutions. Complex heterogeneous media are divided into tetrahedral finite elements inside which the square of the slowness has a simple linear distribution. Arbitrarily complex media may be studied by judicious choice of the elements. We develop appropriate boundary conditions for interfaces with zeroth or first-order velocity discondinuities. Two-point ray tracing is performed by a Newton method based on paraxial ray theory. With a few iterations the ray trajectories through the source and the observer are calculated with the same overall speed as with ray bending. The paraxial method has the additional advantage that it can be used to identify caustics and to separate travel time branches. We used the previously described method to replace the ray tracing algorithm of the HYPO71 program, probably the most commonly used earthquake location program. The ray tracing routine may be easily modified to adapt it to other earthquake location environments. With the modified HYPO71 program, we relocated the aftershock sequence of the November 23, 1980, Irpinia (Italy) earthquake. The introduction of a low-velocity zone proposed by different authors shifts aftershock locations to the NE. This horizontal displacement of hypocenters is not very sensitive to the thickness of the low-velocity zone. ¿ American Geophysical Union 1988