Mullins equation, the fourth-order non-linear partial differential equation describing grain boundary grooving by surface diffusion, is solved by Laplace transform finite difference method (LTFDM). The surface profiles for groove root slopes ranging between 0 and 4.5 are presented. The numerical results show excellent agreement with solutions obtained from Tritscher and Broadbridge’s analytically solvable model when slope at the groove root is not more than 0.7. The shape of the groove is still correct, however, with less accuracy when increasing the groove slope. Other numerical methods are also investigated, including cubic splines and finite difference methods. Our numerical test run of each method shows that, with the same accuracy of results, LTFDM uses time much less than the others.