The second-order saddlepoint approximation (SPA) has been used for component reliability analysis for higher accuracy than the traditional second-order reliability method. This work extends the second-order SPA to system reliability analysis. The joint distribution of all the component responses is approximated by a multivariate normal distribution. To maintain high accuracy of the approximation, the proposed method employs the second-order SPA to accurately generate the marginal distributions of the component responses; to simplify computations and achieve high efficiency, the proposed method estimates the covariance matrix of the multivariate normal distribution with the first-order approximation to the component responses. Examples demonstrate the high effectiveness of the second-order SPA method for system reliability analysis.