A new computational method has been developed for the analysis of three-dimensional large-amplitude motion of liquids with free surfaces in moving containers. The problem is formulated mathematically as a nonlinear initial-boundary value problem under the assumption of irrotational flow of an inviscid fluid. The free surface is moved during each time step in a Lagrangian manner, and its new position is calculated using a forward-time Taylor series expansion. One of the key features of the present method exists in the algorithm to compute the values of Lagrangian time derivatives in the Taylor series. To evaluate these derivatives, a boundary element method is used. The use of the Taylor series expansion has enabled us to employ a variable time-stepping method. The size of a time increment is determined at each time step so that the remainders of the truncated Taylor series should be equal to or less than a given small error limit. Such a variable time-stepping technique has made a great contribution to numerically stable computations. As a numerical example, a swirl motion of liquid in a cylindrical container subjected to a forced horizontal oscillation has been analyzed.