Abstract The application of the method of differential quadrature (Civan et al, 1984, 1988, 1988) is extended for efficient numerical solution of hydrocarbon reservoir models. This method approximates the partial derivatives by a sum of the discrete function values weighted according to a multi-dimensional polynomial. Hence, the differential quadrature method yields accurate solutions without the inherent efficiencies of the conventional methods such as grid point orientation effects, numerical dispersion and oscillations. Differential quadrature approximations can be designed for any order of accuracy without any difficulty. However, it is shown that quadratures of the order of seven to eleven are satisfactory for the solution of highly nonlinear and strongly coupled reservoir model equations. First, the details and a comprehensive theory of the method of differential quadrature are presented. The significant advantages of this method over the conventional finite difference and element methods are demonstrated using model problems having analytic solutions. For this purpose, single phase radial flow problems and the Buckley-Leverett problem are considered for homogeneous and isotropic porous media and constant fluid properties. Second, the method is applied to more realistic cases including heterogeneous and anisotropic media, variable fluid properties and multiphase flow. It is shown that the method of differential quadrature is a rapid, practical and accurate method which circumvents the commonly known difficulties and deficiencies of the conventional methods.