  We extend the branching process based numerical algorithm of Bouchard et al.
[3], that is dedicated to semilinear PDEs (or BSDEs) with Lipschitz
nonlinearity, to the case where the nonlinearity involves the gradient of the
solution. As in [3], this requires a localization procedure that uses a priori
estimates on the true solution, so as to ensure the well-posedness of the
involved Picard iteration scheme, and the global convergence of the algorithm.
When, the nonlinearity depends on the gradient, the later needs to be
controlled as well. This is done by using a face-lifting procedure. Convergence
of our algorithm is proved without any limitation on the time horizon. We also
provide numerical simulations to illustrate the performance of the algorithm.
