This paper summarizes the methods that have been used to develop fully implicit and highly implicit reservoir simulators. The techniques developed emphasize simplicity and clarity. The governing equations, written in fully implicit form are treated as functions. These coupled nonlinear functions are solved by Newtonian iteration, the derivatives of the Jacobian matrix being evaluated numerically. A procedure is described for the efficient evaluation of the Jacobian. A technique is presented for avoiding variable substitution in variable bubble point problems. This is accomplished by means of a novel pseudo solution gas formulation. An efficient method for including two-point upstream into fully implicit simulators, without increasing the Jacobian bandwidth, is described. The scheme is based on a fully implicit upstream mobility plus an explicit correction term derived from the plus an explicit correction term derived from the two-point mobility formulation. In addition a new method called centralized upstream is given. The new method has less truncation error then two-point upstream, and is little more difficult to implement. A modified Crank-Nicholson method is described for reducing the time truncation error associated with large time steps in implicit models. The modified version does not exhibit the usual stability problems associated with the Crank-Nicholson method on nonlinear problems. A section on well models shows how a simple enhancement of conventional well models can be used to account for the position of a well within a grid It is also shown how to introduce a multiblock completion into an implicit simulator without troying the block-banded form of the Jacobian, or without increasing the bandwidth. An efficient iterative solution method is described for the sparse block-structured Jacobians. The method extends the SIP factorization to block form and uses minimizations and orthogonalizations to accelerate the convergence rate. Several applications and examples are provided. A water coning example is included to show the stability and convergence characteristics of the Newtonian iteration. A vertical cross-section oil resaturation problem with variable bubble point is used to illustrate the pseudo solution gas formulation. The modified Crank-Nicholson method is applied to a one-dimensional Buckley-Leverett type problem. problem Introduction It is well known that fully implicit methods are remarkably stable and can tolerate much larger time steps than those used in explicit formulations. However, fully implicit methods are not widely applied in large scale numerical problems such as oil reservoir simulation. The larger matrix bandwidth associated with fully implicit approach demands more computation effort per time step and computer core memory compared with explicit methods. Other considerations such as time truncation error associated with the large time step size and the difficulties in the implementation of higher order methods to reduce spatial truncation have also made the fully implicit method less attractive. A collection of techniques for resolving these difficulties in fully implicit simulation is contained in this paper. GENERAL CHARACTERISTICS OF NEWTONIAN ITERATION The Newton algorithm is a well known iterative process for the numerical solution of equations. process for the numerical solution of equations.