This paper presents a novel reduced-order model (ROM) of thrust for an efficiently flapping airfoil using system identification method. A NACA0012 airfoil pitching and plunging at a low Reynolds number of 40,000 is used to test the ROM. Unlike conventional aerodynamic models which introduce the airfoil displacements directly as inputs, this study utilizes the quadratic terms of displacements as inputs to overcome the frequency-doubling effect of propulsion forces over the oscillation frequency. The autoregressive with exogenous input (ARX) model is adopted to construct mappings between the input and output data. Meanwhile, a heuristic searching strategy is applied for sensitivity analysis of the input variables and the optimal input-vector is determined. The ROMs are then validated in the time domain by comparing their predicted thrust responses with those of CFD simulations under either harmonic or random excitations. Results show that the proposed ROMs can accurately predict the thrust responses of a flapping airfoil with arbitrary motions from moderate to small oscillation amplitudes where a leading edge vortex does not develop, while the computational cost can be reduced by nearly 2 orders of magnitude compared to the high-fidelity CFD simulation method. Finally, the validity of ROMs is mostly clearly shown by using them for propulsive characteristic analysis of a flapping airfoil. Excellent qualities of the ROMs indicate that they can be used for flapping mode optimization and flapping flight control in future research.