Abstract <p>In the present paper, we propose a numerical method for modeling the downstream propagation of optimal disturbances in compressible boundary layers over three-dimensional aerodynamic configurations. At each integration step, the method projects the numerical solution of governing equations onto an invariant subspace of physically relevant eigenmodes; and the numerical integration is performed along the lines of disturbance propagation. The propagation of optimal disturbances is studied in a wide range of parameters for two configurations: a boundary layer over a swept wing of finite span, and a boundary layer over a prolate spheroid. It is found that the dependence of the disturbance energy amplification on the spanwise wavenumber has two local maxima. It is discussed how to combine the developed method with the modern approaches, which are designed to predict the onset of laminar–turbulent transition using the e^<i>N</i>-method.</p>