Parametric partial differential equations are commonly used to model physical systems. They also arise when Wiener chaos expansions are used as an alternative to Monte Carlo when solving stochastic elliptic problems. This paper considers a model class of second order, linear, parametric, elliptic PDE's in a bounded domain D with coefficients depending on possibly countably many parameters. It shows that the dependence of the solution on the parameters in the diffusion coefficient is analytically smooth. This analyticity is then exploited to prove that under very weak assumptions on the diffusion coefficients, the entire family of solutions to such equations can be simultaneously approximated by multivariate polynomials (in the parameters) with coefficients taking values in the Hilbert space [Formula: see text] of weak solutions of the elliptic problem with a controlled number of terms N. The convergence rate in terms of N does not depend on the number of parameters in V which may be countable, therefore breaking the curse of dimensionality. The discretization of the coefficients from a family of continuous, piecewise linear finite element functions in D is shown to yield finite dimensional approximations whose convergence rate in terms of the overall number N dof of degrees of freedom is the minimum of the convergence rates afforded by the best N-term sequence approximations in the parameter space and the rate of finite element approximations in D for a single instance of the parametric problem.