A new analytic element solution has been derived for steady two-dimensional groundwater flow through an aquifer that contains an arbitrary number of elliptical inhomogeneities. The hydraulic conductivity of each inhomogeneity is homogeneous and differs from the conductivity of the homogeneous background. In addition to elliptical inhomogeneities, other elements (such as wells and line sinks) may be present. The method is based on a separable form of the solution for Laplace's differential equation in elliptical coordinates. The piezometric head and the stream function, expressed as continuous spatial functions (as components of the complex potential), may be obtained up to machine accuracy regardless of the shape, size, orientation, and conductivity of the elliptical inhomogeneities. Components of the discharge vector are expressed in a similar manner, using a complex discharge function. Problems with 10,000 or more inhomogeneities can be solved using parallel computing on distributed memory supercomputer clusters. Two examples are included to demonstrate the precision and capabilities of the method. The second example is used to perform a preliminary study of contaminant transport in a highly heterogeneous formation of lognormal conductivity distribution. The results of the transport study are compared with recent theoretical and numerical results that are based on circular inhomogeneities. The new results with elliptical inhomogeneities confirm the findings based on circular inhomogeneities, including a long dispersion setting time and a zero value of the asymptotic transverse dispersivity.