We propose a global convergent numerical method to reconstruct the initial condition of a nonlinear parabolic equation from the measurement of both Dirichlet and Neumann data on the boundary of a bounded domain. The first step in our method is to derive, from the nonlinear governing parabolic equation, a nonlinear systems of elliptic partial differential equations (PDEs) whose solution yields directly the solution of the inverse source problem. We then establish a contraction mapping-like iterative scheme to solve this system. The convergence of this iterative scheme is rigorously proved by employing a Carleman estimate and the argument in the proof of the traditional contraction mapping principle. This convergence is fast in both theoretical and numerical senses. Moreover, our method, unlike the methods based on optimization, does not require a good initial guess of the true solution. Numerical examples are presented to verify these results.