A new stable and accurate three-dimensional (3-D) finite-element scheme based on the split-operator technique, is described for the scalar propagation analysis of optical waveguides with arbitrary cross section and inhomogeneous nonlinear dielectrics, including loss or gain. Recently, serious instability and inaccuracy problems were pointed out for the conventional 3-D split-operator finite-element scheme. The present algorithm overcomes these problems by preserving all the attractive features of the split-operator technique. Due to the finite-element discretization, computer resources can be optimized by exploiting the use of adaptive meshes and sparse-matrix techniques.