A new method is presented to compute gravity anomalies that result from 2-D bodies with variable density contrast. The cross-section of a body is approximated by a polygon. Density is assumed to vary as any order of polynomial function with depth and lateral position. Results calculated by the proposed method for models with variable density contrast compare well with other methods. Liaohe basin, northeast China, is modeled from field gravity anomalies using the formulas given, showing the method is valid and effective.