When we use penalty method to solve the pure Neumann problem, the associated stiffness matrix we get is ill-conditioned, which leads to the unstable computation. In this paper, we design an iterative penalty method for the problem and prove some error estimates. In our algorithm, we can use a not very small penalty parameter to avoid the unstable computation. Numerical examples are given to show the algorithm is very effective and powerful. (c) 2006 Elsevier Inc. All rights reserved.