We develop an ELLAM-MFEM solution procedure for the numerical simulation of compressible fluid flows in porous media with point sources and sinks. An Eulerian–Lagrangian localized adjoint method (ELLAM), which was previously shown to outperform many widely used and well-regarded methods in the context of linear transport partial differential equations, is presented to solve the transport equation for concentration. Since accurate fluid velocities are crucial in numerical simulations, a mixed finite element method (MFEM) is used to simultaneously solve the pressure equation as a system of first-order partial differential equations for the pressure and mass flow rate. This minimizes the numerical difficulties occurring in standard methods caused by differentiation of the pressure and then multiplication by rough coefficients. Computational experiments show that the ELLAM-MFEM solution procedure can accurately simulate compressible fluid flows in porous media with coarse spatial grids and very large time steps, which are one or two orders of magnitude larger than those used in many numerical methods. The ELLAM-MFEM solution technique symmetrizes the governing partial differential equations, and greatly reduces or eliminates non-physical oscillation and/or excessive numerical dispersion present in many large-scale simulators that are widely used in industrial applications. It conserves mass and treats boundary conditions in a natural manner. It can treat large adverse mobility ratios, discontinuous permeabilities and porosities, anisotropic dispersion in tensor form, compressible fluid, heterogeneous media, and point sources and sinks.