Smoothing thin plate splines, a fitting technique based on a rigorous roughness penalty approach, have been recently investigated as a promising tool for bivariate interpolation of aerodynamic data. In this paper, this technique is implemented and extended to multivariate fitting. In particular, the method is applied for estimating the aerodynamic polars of well-known two-dimensional symmetrical and nonsymmetrical airfoils as functions of some geometric parameters describing the airfoil shape and a further variable defining the flow regime (either the Mach or the Reynolds number). Therefore, the simultaneous influence of five independent variables on three responses (lift, drag, and pitching moment coefficients) is investigated. To this purpose, a large database is generated via numerical simulations (using a validated flow solver) containing all information required to build a reliable response surface. Then, the model is built and its performance validated by performing queries on complete aerodynamic polars at various flow regime conditions of a series of airfoils not included into the database. Results show a very good matching between predicted and calculated curves, thus demonstrating the remarkable predictive capability of the implemented tool.