We propose a method that significantly improves the accuracy of the level set method and could be of value for numerical solutions of differential equations in general. Level set methods use a level set function, usually an approximate signed distance function, Φ, to represent the interface as the zero set of Φ. When Φ is advanced to the next time level by an advection equation, its new zero level set will represent the new interface position. But the non-zero curvature of the interface will result in uneven gradients of the level set function which induces extra numerical error. Instead of attempting to reduce this error directly, we update the level set function Φ forward in time and then backward to get another copy of the level set function, say Φ 1 . Φ 1 and Φ should have been equal if there were no numerical error. Therefore Φ−Φ 1 provides us the information of error induced by uneven gradients and this information can be used to compensate Φ before updating Φ forward again in time.