The joint frequency spectrum (JFS) is a fundamental summary for genetic polymorphism of samples from multiple populations, and is useful for population genetic inference of ancient demography. Analytical theory and methodology of the JFS have been developed using coalescent theory, but are only practicable for small samples. When the sample size is large (e.g., n>50), the computation of JAFS becomes numerically intractable. We present accurate approximation for the JFS for large samples from one or multiple populations. The novelties include: (1) the exact formulas of the JFS are in simple form and computationally efficient for large samples without the numerical intractability existing in former studies; (2) arbitrarily time-varying population size can be modeled. The accuracy of the results is demonstrated by comparing to coalescent simulations. The work provides an useful tool for JFS-based demography inference using large-sample genomic sequencing data.