Perturbation theory (PT) calculation of large-scale structure has been used to interpret the observed non-linear statistics of large-scale structure at the quasi-linear regime. In particular, the so-called standard perturbation theory (SPT) provides a basis for the analytical computation of the higher-order quantities of large-scale structure. Here, we present a novel, grid-based algorithm for the SPT calculation, hence named GridSPT, to generate the higher-order density and velocity fields from a given linear power spectrum. Taking advantage of the Fast Fourier Transform, the GridSPT quickly generates the nonlinear density fields at each order, from which we calculate the statistical quantities such as non-linear power spectrum and bispectrum. Comparing the density fields (to fifth order) from GridSPT with those from the full N-body simulations in the configuration space, we find that GridSPT accurately reproduces the N-body result on large scales. The agreement worsens with smaller smoothing radius, particularly for the under-dense regions where we find that 2LPT (second-order Lagrangian perturbation theory) algorithm performs well.