The laminar-turbulent transition of boundary layers is accompanied by an increase in friction drag. Laminar airfoils are a promising technology to reduce the emission of greenhouse gasses caused by the aerospace industry. The transition scenario in airliners is dominated by stationary crossflow vortices whose evolution is strongly nonlinear. This transition is preceded by the eigenmode growth of instabilities, rise of secondary instabilities and development of turbulent spots that eventually lead to a fully turbulent boundary layer. Exponential growth starts at infinitesimal amplitudes at which linear methods, e.g. the Orr-Sommerfeld or Linear Parabolized Stability Equations (LPSE), can be exploited to assess the growth of instabilities. As the energy of these instabilities increases, interactions become important and a nonlinear method capable of calculating wave-triad interactions should instead be used to correctly predict their evolution. This thesis concerns the development of a tool capable of solving the Nonlinear Parabolized Stability Equations (NPSE) and introduces a new method for the introduction of harmonics during the simulation that follow from nonlinear interactions. This method accounts for nonparallel and nonlinear effects in an Inhomogeneous LPSE (ILPSE) framework by modeling the interactions as a forcing term in this equation. The amplitude of the mode prior to its introduction is modeled via the nonlinear estimation of the growth rate. An investigation in mode introduction schemes proved that the introduction of modes severely affects the results from stability analyses. The ILPSE method reduces initial transients by correctly accounting for a harmonic’s history using a growth rate estimate from the other modes in the system. The improved introduction scheme severely increases the accuracy of the introduction amplitude and growth of higher harmonics.