An apparatus and method are provided for solving a non-linear S-shaped function F=f(S) which is representative of a property S in a physical system, such saturation in a reservoir simulation. A Newton iteration (T) is performed on the function f(S) at S.sup.v to determine a next iterative value S.sup.v+1. It is then determined whether S.sup.v+1 is located on the opposite side of the inflection point S.sup.c from S.sup.v. If S.sup.v+1 is located on the opposite side of the inflection point from S.sup.v, then S.sup.v+1 is set to S.sup.l, a modified new estimate. The modified new estimate, S.sup.l, is preferably set to either the inflection point, S.sup.c, or to an average value between S.sup.v and S.sup.v+1, i.e., S.sup.l=0.5(S.sup.v+S.sup.v+1). The above steps are repeated until S.sup.v+1 is within the predetermined convergence criteria. Also, solution algorithms are described for two-phase and three-phase flow with gravity and capillary pressure.