An improved matrix method for the analysis and design of nonlinear directional couplers (NLDC) with Kerr-like coupling media is presented. The method represents an extension of the original matrix method that has been used to study linear waveguides with inhomogeneous refractive index profiles and three-layer nonlinear waveguides. The method is based on an iterative averaging algorithm that calculates the average values of the dielectric constant and the field amplitude in each stratified nonlinear layer. This numerical method is applied to compute the power-dependent characteristics of symmetrical and nonsymmetrical multiple quantum-well (MQW) NLDCs. Saturation of the refractive index in the MQW coupling medium is taken into account and, for a symmetrical structure, the numerical results are shown to be in agreement with published experimental data.