A two-dimensional computational model has been developed to study the evolution and breakup of a viscous laminar liquid jet, using a boundary-fitted curvilinear coordinate system. A system of elliptic partial differential equations for coordinate transformations has been developed to map the moving boundaries’ physical domain of the jet to a simple rectilinear computational domain. The equations developed for the model comprise the transformed two-dimensional unsteady Navier–Stokes equations for the liquid jet, grid velocity equations, kinematic boundary conditions, and the Geometric Conservation Law. The resulting systems of equations are solved using an implicit finite difference scheme. Effects of inflow oscillation magnitude, wave number, Weber number, and Reynolds number on the breakup process of jets have been studied. The model predicts the instantaneous shape of the jet surface, formation of the main and satellite drops, and the breakup length and time. These results are compared with available experimental data. The comparisons show a good agreement between measured and computed values of drop sizes and breakup lengths for different Reynolds and Weber numbers. However, at a relatively high Reynolds number of 1,254, the model slightly overpredicts the main drop sizes and underpredicts the satellite drop sizes at a wave number of 0.4. At a low Reynolds number of 587, the model overpredicts the main drop sizes at a lower wave number of 0.3. Moreover, the model underpredicts the satellite drop sizes at a lower wave number of about 0.4 and overpredicts the satellite drop sizes at a wave number of 0.8.