The design of natural-laminar-flow airfoils is demonstrated by high-fidelity, multipoint, aerodynamic shape optimization capable of efficiently incorporating and exploiting laminarturbulent transition. First, a two-dimensional Reynolds-averaged Navier-Stokes (RANS) flow solver has been extended to incorporate an iterative laminar-turbulent transition prediction methodology. The natural transition locations due to Tollmien-Schlichting instabilities are predicted using the simplified e N envelope method of Drela and Giles or alternatively, the compressible form of the Arnal-Habiballah-Delcourt criterion. The boundarylayer properties are obtained directly from the Navier-Stokes flow solution, and the transition to turbulent flow is modeled using an intermittency function in conjunction with the Spalart-Allmaras turbulence model. The RANS solver is subsequently employed in a gradient-based sequential quadratic programming shape optimization framework. The laminar-turbulent transition criteria are tightly coupled into the objective and gradient evaluations. The gradients are obtained using a new augmented discrete-adjoint formulation for non-local transition criteria. The aerodynamic design requirements are cast into a multipoint design optimization problem. A composite objective is defined using a weighted integral of the operating points. A Pareto front is also formed to study and quantify off-design performance. The proposed framework is applied to the single and multipoint optimization of subsonic and transonic airfoils, leading to robust natural-laminar-flow designs.