High fidelity models for uncertainty quantification of large structures in finite element framework are computationally exhaustive. Thus, there is a constant demand for efficient algorithm that uses optimal computational cost without compromising with the quality of the end results. With this in view, present study aims to develop a sequentially evolving stochastic response surface using Hermite family of orthogonal polynomials whose support points are generated in sparse grid framework. Using the values of the original model at these support points, unknown coefficients of the stochastic response surface are optimized by moving least square technique. It helps to reduce the number of original function evaluations to determine the coefficients as compared to other deterministic or random sampling techniques. Besides sparse grid scheme for support point generation, they are also populated sequentially as the optimization progresses in every iteration. The uniqueness of the proposed scheme is its adaptability by changing the order of the polynomials and the level of the sparse grid to minimize the overall computational cost. Multiple optima present in the original response can be identified by introducing additional penalty functions whenever they are required. Once the global response surface is ready, Monte-Carlo simulation or its advance version (e.g. Latin Hypercube Sampling) is adopted to generate the probability density functions for the output variables. Numerical studies are presented to prove the efficiency and accuracy of the proposed scheme as compared to other techniques available in the literature. © 2017 The Authors.