A numerical scheme based on the Galerkin’s finite element method is used to analyze nonisothermal flows of Bingham plastics. These fluids exhibit yield stresses, below which the material does not flow. This condition is enforced by introducing the bi-viscosity model in which the plug behaviour is approximated by a highly viscous fluid. The bi-viscosity model approaches the ideal Bingham plastic model if the pre-yield viscosity becomes large. This numerical method is applied to solve three problems namely the Graetz-Nusselt, recirculating flow, and sinusoidal channel flow problems for various values of the yield stress. The size of the plug as well as the velocity and temperature fields correspond closely with values available in literature.