A spectral relaxation method used with bivariate Lagrange interpolation is used to find numerical solutions for the unsteady three-dimensional flow problem of an Oldroyd-B fluid with variable thermal conductivity and heat generation. The problem is governed by a set of three highly coupled nonlinear partial differential equations. The method, originally used for solutions of systems of ordinary differential equations is extended to solutions of systems of nonlinear partial differential equations. The modified approach involves seeking solutions that are expressed as bivariate Lagrange interpolating polynomials and applying pseudo-spectral collocation in both independent variables of the governing PDEs. Numerical simulations were carried out to generate results for some of the important flow properties such as the local skin friction and the heat transfer rate. Numerical analysis of the error and convergence properties of the method are also discussed. One of the benefits of the proposed method is that it is computationally fast and gives very accurate results after only a few iterations using very few grid points in the numerical discretization process.