In this paper, a least-squares spectral element method for parabolic initial value problem for two space dimension on parallel computers is presented. The theory is also valid for three dimension. This method gives exponential accuracy in both space and time. The method is based on minimization of residuals in terms of the partial differential equation and initial condition, in different Sobolev norms, and a term which measures the jump in the function and its derivatives across inter-element boundaries in appropriate fractional Sobolev norms. Rigorous error estimates for this method are given. Some specific numerical examples are solved to show the efficiency of this method. © 2017, Springer Science+Business Media New York.