||Drop formation is of great importance in scientific research as well as in industrial applications. Many researches have been carried out experimentally or theoretically. Numerical simulation serves as a good complementary of the experimental investigation and theoretical analysis. In this thesis, we numerically simulate the dynamics of drop formation in two-phase flows. We first investigate drop formation dynamics in a co-flowing fluid with different density and viscosity based on phase field model. We simulate a three-dimensional (3D), axi-symmetric system consisting of Cahn-Hilliard equation and Navier-Stokes equation under cylindrical coordinate using finite difference method. A new fractional time-stepping technique proposed recently is applied to solve the pressure. Only a Poisson equation needs to be solved per time step, just like what one needs to do for two-phase fluids with equal density. The Cahn-Hilliard equation is solved implicitly using convex-splitting method. It is shown that the limiting length, the volume of the drop, as well as the shape of the drop at the breakup during steady state are affected by the dimensionless parameters. The numerical results are compared with Zhang’s experiment [Z99b], the shape of the bubble during its evolution matches well with the experiment. We test the dripping to jetting transition for 5 different cases, the numerical results are found to be in qualitative agreement with the experiment carried out by Utada et. al [UFNSW07]. We then study drop formation of an incompressible Newtonian fluid flowing through a vertical tube into air. This phenomenon can be described by a one-dimensional (1D) system of partial differential equations (PDE) which develops singularity. This singularity has similarity solutions that can be obtained by solving a 1D system of ordinary differential equations (ODE). We solve the 1D PDE system adaptively using iterative grid distribution method. High resolution near the singularity region is achieved. The solutions of the 1D ODE system are solved numerically both before and after the breakup, based on the successful derivation of general expressions for all the coefficients of the similarity functions. It is found that, the numerical solutions of the PDE system collapse onto the numerical solutions of the ODE system as the time distance of the singularity approaches to zero.