This paper describes the application of Newton's method to low Mach number steady and unsteady laminar diffusion flames, which are characterized by low flow speed and highly variable density. The newly-emerged vorticity-velocity formulation of the Navier-Stokes equations is used for both steady and unsteady compressible flows to avoid staggered mesh discretization. The nonlinear Navier-Stokes equations are discretized using finite difference method, and a secondorder backward Euler scheme is applied for the time derivatives. Central difference is used for diffusion terms to achieve better accuracy, and a monotonicity-preserving upwind difference is used for convective ones. We use an unequal-sized single grid mesh for unsteady flow and a three level multigrid method for steady flow. The coupled nonlinear system is solved via the damped Newton's method for both steady and unsteady flows. The Newton Jacobian matrix is formed numerically, and the resulting linear system is ill-conditioned and is solved by the iterative solver Bi-CGSTAB with a Gauss-Seidel preconditioner.