In this research, the plunging motion of an airfoil by a numerical method based on finite volume in different Reynolds numbers is simulated and the thickness effect, amplitude and reduced frequency on the aerodynamic coefficients are investigated. In this process, SIMPLEC algorithm, implicit solver, high order scheme and dynamic mesh method is used in unsteady simulation and the flow is supposed viscous, incompressible and laminar. Simulations are in three Reynolds 1000, 11000 and 50000, respectively, in accordance with the flight of the insects, small birds and pigeons are done in two amplitudes and three reduced frequencies. The simulation results are compared with published data to confirm the validity of research. This comparison shows comprisable agreement. Pressure distribution and Vortex shedding around airfoils show that the thickness of the airfoils delays vortex shedding and changes time-averaged thrust coefficient. Reduced frequency and amplitude of oscillation are two important parameters in this simulation, but the reduce frequency is more effective than amplitude. The response surface methodology (RSM) was used to optimize the plunging airfoil. Optimization shows that airfoil with 0.29% thickness, 3.08 reduced frequency and 0.5 dimensionless oscillating amplitude produce maximum trust coefficient.