Abstract: A high-order spectral difference (SD) scheme based cost effective optimization approach is proposed to optimize the kinematics of flapping airfoils for maximum propulsive efficiency. Specifically, a gradient based optimization algorithm is coupled with a high-order SD Navier-Stokes solver to search for the optimal kinematics of a series of NACA 4-digit airfoils undergoing combined plunge and pitch motion on a coarse mesh. Then the high-order SD solver is used to capture the detailed vortex structures associated with the optimal kinematics of the flapping flight on a fine mesh. The proposed numerical optimization framework is used to study the effects of airfoil thickness, Reynolds number, and pitching pivot location for optimal cruise flight. The flow physics behind optimal airfoil oscillation is explained by examining the corresponding flow fields, aerodynamic forces and the variation of effective angle of attack.