To investigate the hydrodynamic performance of flexible flapping foils and provide guidance for the design of bio-inspired propulsors, numerical simulations are used to investigate the unsteady flow around a three-dimensional (3D) flexible foil based on an in-house developed immersed boundary method CFD code. The effect of complex moving deformation boundary on the flow on a fixed Cartesian gird is imposed by a boundary condition reconstruction algorithm based on ghost-cells. A numerical method calculating propulsive efficiency which can apply to both rigid and flexible flapping foils is proposed in views of energy. Moreover, the effect of various parameters on the propulsive performance of the 3D flexible foil has been systematically investigated. The results indicate that the propulsive efficiency of the flexible flapping foil is higher than its rigid counterpart when deformation parameter ε lies between 2.0 and 3.5 together with chordwise deformation coefficient δ of 0.1. The variation of the flexible foil performance with deformation parameter depends on the vortex dynamics underlying the force production, including the control of the trailing-edge deformation on the vortex shedding and the mechanism of spanwise transport of vorticity. The mean thrust of the flexible flapping foil increases monotonically with the dimensionless frequency k, which is consistent with the experimental phenomenon that fish propelled by pectoral fins enhance the thrust through the increase of the flapping frequency of pectoral fins. There exists an optimal k to get maximum propulsive efficiency of the flexible foil, and the simulation reproduces the characteristic of fish propelled by caudal fin in experiments, i.e., they cruise at a Strouhal number tuned for optimal efficiency.