Polymer mold filling is a non-Newtonian, non-isothermal, non-steady process with free surface. It is well known that the fountain flow effect in the flow front region and the solidified layer near the cold mold wall have significant effect on the mold filling process and thus the mechanical properties of the molded parts. In this study, a non-isothermal, non-steady mold filling process including the solidified layer as well as the fountain flow effect was analyzed numerically using the finite element method. Shear rate dependency in the viscosity was taken into account by employing the Cross model. Single-step algorithm was used for the numerical integration in time and the flow front was advanced by the explicit scheme. The resulting system of algebraic equations was solved using the Picard iteration technique. In order to verify the validity of the numerical scheme, energy balance and mass conservation were checked. It was found that the energy and mass balances were satisfied within reasonable accuracies. The importance of the non-steady term was assessed by analyses done in a similar manner by assuming that the process is quasi-steady. Two results thus obtained were compared. The quasi-steady assumption was found to be insufficient in describing the inherent non-steady effect during mold filling.