This article presents a p-version least-squares finite-element formulation for axisymmetric incompressible Newtonian and non-Newtonian fluid flow with heat transfer. The dimensionless form of the equations describing the fluid motion and heat transfer are cast into a set of first-order coupled partial differential equations involving non-Newtonian stresses and heat fluxes as auxiliary variables. The pressure, velocities, stresses, and heat flares are interpolated using equal-order, C-0 continuity, p-version hierarchical approximation functions. The least-squares error functional is constructed using the system of coupled first-order nonlinear partial differential equations without linearization, approximations, or assumptions. The minimization of this least-squares error functional results in a solution vector {delta} that makes the partial derivatives of the error functional with respect to {delta} a null vector. This is accomplished by using Newton's method with a line search. The article presents implementation of the power law model for non-Newtonian viscosity. The fluid properties are also considered to be a function of temperature. Numerical examples are presented to demonstrate the convergence characteristics and the accuracy of the proposed formulation.