Heat transfer of a power-law non-Newtonian incompressible fluid in channels with porous walls has not been carefully studied using a proper numerical method despite a few constructions of approximate analytic solutions through the similarity transformation and perturbation method for Newtonian fluids (i.e. power-law index being one). In this paper, we propose a finite element method for the thermal incompressible flow equations. The incompressible condition is treated by a penalty formulation. Numerical solutions are validated by comparing them with an approximate analytic solution of the Navier-Stokes equation in the Newtonian fluid case. Then, the method is used to simulate the heat transfer of various power-law fluids. Additionally, unlike previous studies, we allow the thermal diffusivity to be a function of temperature gradient. The effect of different values of the parameters on the temperature and velocity is also discussed in this paper.