We employ an accurate numerical algorithm to simulate two model flames—an unconfined lifted and a confined, coflowing, methane-air jet diffusion flame using detailed chemistry and complex transport models. The algorithm employs Newton's method to obtain the primitive variable solution of the large system of strongly coupled elliptic governing equations. The Newton equations are solved by a block-line tridiagonal method. We employ a global grid refinement technique which equidistributes meshes according to the gradients and curvatures of the solution obtained on the previous mesh and bounds the ratio of adjacent grid step size. The algorithm can be applied to problems ranging from non-reacting flows to reacting flows in two- or three-dimensional configurations. In the unconfined case, the lifted flame and the “triple flame” are both predicted in the numerical solution. The computed solutions agree well with the experimental results. The comparison of the present solutions with the previously reported stream function-vorticity computations shows that the primitive variable numerical approach is robust and yields a more accurate solution due to its implementation of the flow boundary conditions.