The authors propose and analyze a gradient recovery technique for finite element solutions which provides new gradient approximations with high order of accuracy. They modify the method of O. C. Zienkiewicz and J. Z. Zhu [Comput. Methods Appl. Mech. Eng. 101, No. 1–3, 207–224 (1992; Zbl 0779.73078)] by applying a global least-squares fitting to the gradient of finite element approximation and provide a theoretical analysis for the modified Zienkiewicz-Zhu method by establishing a superconvergence estimate for the recovered gradient/flux on general quasi-uniform meshes. They also give numerical results to show that the recovery technique is robust and efficient.