Many problems in power systems depend on parameters, which could be stochastic variables or deterministic system control variables practically, e.g., generation outputs, nodal voltages, etc. Due to the nonlinearity of power systems, the analytical relation between system states and parameters cannot be obtained directly. Using the sampling method to evaluate the influence of parameters on system states is very powerful but time-consuming. One feasible approach is to use polynomial approximations, where the system states are approximately expressed in the form of polynomials in terms of parameters. Galerkin method can be used to identify the approximate solution with high accuracy by solving high-dimensional equations. However, if a large number of parameters are involved, solving these high-dimensional equations becomes a serious challenge. This paper proposes an innovative method for resolving these high-dimensional equations in power systems, which constructs a sequence of decoupled equations to determine the polynomial expansion coefficients. This new approach can provide a local approximation in the form of Taylor expansion at a given operation point. Although the method is general, for simplicity and good readability, we introduce the detailed process in its application to load flow problems. Case studies from 6-, 118-, and 1648-bus system show that the proposed method provides approximation more efficiently than traditional Galerkin method does, and 3-order polynomials can give very accurate results.