This paper presents 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. This paper also attempts to enhance the computational efficiency by presenting an interpolation method used to evaluate matrix elements, which are ascribed to a discretization procedure. In addition to the steady case, numerical examples consist of cases involving a submerged prolate spheroid, which is originally idle and then suddenly moves with a constant speed and a constant acceleration. Also systematically studied is the variation of hydrodynamic forces acting on the spheroid for various Froude numbers and submergence depths.