The major drawbacks of the Hele-Shaw approximation, commonly used today as a means of simplifying the simulation of the injection molding process, are the inherent loss of the ability to predict important physical phenomena and the ambiguity involved in the definition of a midplane. With the objective of eliminating many of these problems, a fully three-dimensional mold filling simulation program has been developed. The governing equations are in terms of Navier-Stokes problems for the viscous, incompressible, nonisothermal, and non-Newtonian fluid. To avoid the simultaneous determination of the coupled velocity and pressure, this article introduces an iterative method that at any given time step solves the components independently. This method reduces memory needs in simulation and enhances the stability of numerical scheme. In addition, this article also presents a mixed implicit and ''up-wind'' scheme to discrete the energy equation. It can overcome the spatial oscillations of temperature in numerical simulation. Good agreements with both analytic solution and the actual processes are found in the current investigation. This method can successfully predict the flow features in injection molding.