Based on the finite difference scheme of Navier-Stokes equations, the dynamic response of underwater bubbles is numerically simulated. A high accuracy numerical method is realized to define the free surface boundaries with severe deformation. An innovative outer boundary condition is proposed to improve the numerical efficiency, which was low in the past applications. First, the flow pattern far from the bubble is assumed to be generated by point-source/point-sink and therefore the pressure can be calculated. If this pressure is used as the boundary condition, the infinite domain problem can be converted to small region simulation. Invalid outer boundary pressure phenomena and numerical oscillation phenomena are inhibited by distinguishing the pressure’s rationality. The method is proved to be of high efficiency and accuracy. Good agreement is obtained between numerical results, Rayleigh-Plesset equation and experiments. Finally, the dynamic response of a bubble in the vicinity of a rigid wall is simulated and the process of jet formation is analyzed. The method can be used for the simulation of the underwater bubble pulsation as well as the prediction of jet loads.