We present a numerical method to evaluate the hydrodynamic forces of translating bodies under a free surface. Both steady and unsteady problems are considered. Analytical and numerical studies are carried out based on the Havelock wave-source function and the integral equation method. Two main problems arising inherently in the proposed solution method are overcome in order to facilitate the numerical implementation. The first lies in evaluating the Havelock function, which involves integrals with highly oscillatory kernels. Particular integration contours leading to non-oscillatory integrands are derived a priori so that the integrals can be evaluated efficiently. The second problem lies in evaluating singular kernels in the boundary integral equation. The corresponding non-singular formulation is derived using some theorems of potential theory, including the Gauss flux theorem and the property related to the equipotential body. The subsequent formulation is amenable to the solution by directly using the standard quadrature formulas without taking another special treatment