A mode vector parabolic equation (MVPE) is derived for the propagation of normal modes in the ocean waveguide. This model, which includes mode coupling, is a generalization of the adiabatic mode PE [M. D. Collins, 2269 (1993)]. The main features of this model are: (1) It is based on an initial value differential equation which propagates normal modes rather than a method that involves matching solutions at interfaces. (2) The solution vector whose components are the amplitude of the modes is not derived from a procedure involving a reference wave number and thus each vector component is locally accurate around each mode. (3) It agrees with previous benchmark solutions. (4) The model is tractable in three dimensions. MVPE has three parts: Spatial derivative of the mode amplitudes, eigenvalue, and coupling coefficients matrices. The latter two are precomputed; the coupling matrix is based on McDonald’s formulation in terms of environmental parameters which is summarized in [McDonald etal., 2357 (1994)]. The system of equations can either be solved by matrix inversion or iteration. To illustrate the method, the solution to the problem of propagation of waves in a two-dimensional wedge is presented in detail.