The dynamic equations of motion for a three-dimensional piping system conveying a harmonically pulsating fluid contain time-varying terms attributable to the fluid pulsating in the pipe. In this study, finite element formulation for this three-dimensional piping system was performed. The stiffness and damping matrices in the finite element modeling vary according to time because of the effects of the harmonically pulsating fluid. The frequency-domain method based on eigenvalue analysis cannot be used in this kind of problem. Conventional numerical time-domain methods require substantial computational efforts. An efficient numerical method to predict the steady-state time response of the piping system was presented. In this method, simultaneous equations were constructed by comparing the coefficients of a Taylor series expansion instead of directly solving the problem in the time domain. The accuracy and efficiency of this method were validated by comparison with a conventional numerical integration method.