This paper discusses a direct method for numerical computation of the local quadratic approximations of (n - 1)-dimensional stable manifolds for the unstable equilibria on the stability boundary in large systems. The algorithm proposed here does not require time-consuming normal form or symbolic computations. The numerical procedure is illustrated on transient stability problems in simple power system models.