A new numerical simulation code is developed for elucidating the behavior of a bubble being based on the Navier-Stokes equation. The capturing of the interface between liquid and gas is implemented by the modified density function method significantly suppressing numerical diffusion in the rectangular coordinate system. The density function method, which eliminates the efforts of generating grids fitted to the interface, is shown to capture the complex 3-dimensional deformation of a bubble with sufficient degree of accuracy. Computational results of a rising bubble at high Reynolds number greater than 1000 showed the unsteady motions such as those with zigzag or spiral trajectory and the various deformation of bubble shape like a spherical cap. It is also clarified that the unsteady motion is caused by the asymmetric formulation of longitudinal vortex behind the bubble. This code can be applied to the case of bubbles in a boundary layer to investigate the effects of bubbles on the boundary layer properties.