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 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 component responses. Examples demonstrate the high effectiveness of the second order SPA method for system reliability analysis.