Computational methods for additional stresses in foundations are an important topic in soil mechanics. This paper develops a new numerical method for additional elastic fields in layered and non-homogeneous foundations. The techniques adopted primarily involve the use of a two-dimensional integration to integrate the fundamental solutions of a transversely isotropic and layered halfspace under point loads over the distributed pressure. To analyze non-uniform loads over an arbitrary area, the loading area needs to be discretized into a finite number of quadrilateral elements. The regular and singular integrals in the expressions of elastic fields in foundations are accurately calculated. Numerical examples show that the proposed method can obtain highly accurate results and can analyze the foundations with the elastic parameters arbitrarily varying with depth. The proposed method is explicitly and concisely expressed for effective teaching and learning.