In this paper, lift and drag coefficients were numerically investigated using NUMECA software in a set of 4-digit NACA airfoils. Two metamodels based on the evolved group method of data handling (GMDH) type neural networks were then obtained for modeling both lift coefficient (CL) and drag coefficient (CD) with respect to the geometrical design parameters. After using such obtained polynomial neural networks, modified non-dominated sorting genetic algorithm (NSGAII) was used for Pareto based optimization of 4-digit NACA airfoils considering two conflicting objectives such as (CL) and (CD). Further evaluations of the design points in the obtained Pareto fronts using the NUMECA software showed the effectiveness of such an approach. Moreover, it was shown that some interesting and important relationships as the useful optimal design principles involved in the performance of the airfoils can be discovered by the Pareto-based multi-objective optimization of the obtained polynomial meta-models. Such important optimal principles would not have been obtained without using the approach presented in this paper.