An effective method for computing the large scale time-delay stability region of a power system is presented. Firstly, a transformation is introduced to the original time-delay equations. Secondly, an iterative process in a limited subspace is performed to find the core boundary of stability region, which is then cloned, panned and converted into the time-delay space to locate the boundary of power system stability region. Thirdly, a four-point interpolation process is used to distinguish the stable and unstable regions. Finally, some typical systems are employed to validate the presented method. It is revealed that the method can compute the stability region of power systems with time-delay in a large time-scale precisely.