A number of flow cases of practical significance exhibit a predominant spatial direction, along which the mean properties of the flow field vary slowly while having fast variations on the cross-sectional planes. This property is taken into account when the three- dimensional parabolized stability equations (PSE-3D) are derived. These equations represent the most efficient approach for the solution of the instability problem of such flows. In this work, the linear PSE-3D are extended to predict the nonlinear development of perturbations in this kind of complex three-dimensional flows. The newly developed method is formulated and verified for different flow problems of interest. Firstly, it has been verified by computing the evolution of linear and nonlinear Tollmien- Schlichting waves in Blasius boundary layer, showing excellent agreement with traditional nonlinear PSE predictions. Also, the evolution of optimal streaks is simulated and compared against direct numerical simulations. Finally, the nonlinear development of stationary crossflow instabilities in a three-dimensional boundary layer is monitored using a non-orthogonal coordinate system to follow the instability trajectory, showing again a very good agreement with PSE results.