Phase change problems, occurring e.g. in melting, casting and freezing processes, are often characterized by a very narrow transition zone with very lareg changes in heat capacity and conductivity. This leads to problems in numerical procedures, where the transition zone propagates through a mesh fixed in space. An iterative time-stepping method has been developed within a finite element context which improves upon present methods in three ways: a flux based evaluation of the conductivity matrix replaces the traditional temperature based method, a predictor is introduced for the properties in the mid-point of the current time step and full energy balance is restored in each iteration by adjusting the local temperature. Examples demonstrate improved accuracy, a reduction of the numerical disturbances introduced by the passage of the transition zone and a smaller number of timesteps needed compared with other methods currently in use.