Recent studies indicate that flapping-foil biomimetic systems can significantly improve ship propulsion in waves. In most cases horizontal wing(s) are used performing vertical oscillatory motion induced by ship heaving and pitching responses in waves, while foil's rotational motion is actively controlled. Aiming to the detailed investigation of the free-surface effects, a GPU-accelerated Boundary Element Method (GPU-BEM) has been developed for the hydrodynamic analysis of 3D lifting bodies in non-linear waves with an angle of incidence and motions. The nonlinear free-surface conditions are satisfied at the exact location of the free-surface boundary. A nonlinear pressure-type Kutta condition is applied at the trailing edge and the dynamic evolution of the free vortex wake is obtained by a time-stepping method. The free-surface boundary conditions and Kutta condition are used to set up the nonlinear and nonlocal evolution equations. For the treatment of the dynamical system, the present BEM is coupled to the dynamical equations and it is integrated in time domain with an iterative scheme. Results are presented for a few subproblems of special interest, such as the wave resistance and the enforced radiation problems of smooth immersed bodies and compared with other methods for verification. Subsequently, the unsteady biomimetic thruster is examined in unbounded domain, beneath the free surface in the absence of incident waves, and in the presence of obliquely incident waves, illustrating the performance of the system. The present method can be exploited for the systematic investigation and the design of flapping-foil thrusters for augmenting ship propulsion in waves.