A numerical optimization algorithm based on the steepest decent along the variation of the optimization function is implemented for maximizing the thrust and/or propulsive efficiency of a single flapping airfoil. Unsteady, low speed laminar flows are computed using a Navier-Stokes solver on moving overset grids. The flapping motion of the airfoil is described by a combined sinusoidal plunge and pitching motion. Optimization parameters are taken to be the amplitudes of the plunge and pitching motions, and the phase shift between them. Computations are performed in parallel in a work station cluster. The numerical simulations show that high thrust values may be obtained at the expense of reduced efficiency. For high efficiency in thrust generation, the induced angle of attack of the airfoil is reduced and large scale vortex formations at the leading edge are prevented.