Due to the roaring power dissipation and gaining popularity of 3D integration, thermal dissipation has been a critical concern of modern VLSI design. The availability for chip-level 3D thermal analysis is crucial for thermal integrity analysis. In this paper, we formulated an analytical Green function solution of the 3D temperature distribution in a multi-layer integrated circuit substrate. The proposed analytical solution differs from the previously reported results in two ways: (a) This is a 3D temperature distribution solution, while only 2D solutions are reported in the past; and (b) a rectangular coordinate system is used rather than the cylindrical coordinates used in the past which leads to a more proper and accurate representation to the VLSI geometry. Experimental results demonstrate the speed advantage of our approach and the accuracy within the error 0.5% when compared with commercial computational fluid dynamics software ANSYS[13].