  We present an iterative root finding method for harmonic mappings in the
complex plane, which is a generalization of Newton's method for analytic
functions. The complex formulation of the method allows an analysis in a
complex variables spirit. For zeros close to poles of $f = h + \bar{g}$ we
construct initial points for which the harmonic Newton iteration is guaranteed
to converge. Moreover, we study the number of solutions of $f(z) = \eta$ close
to the critical set of $f$ for certain $\eta \in \mathbb{C}$. We provide a
Matlab implementation of the method, and illustrate our results with several
examples and numerical experiments, including phase plots and plots of the
basins of attraction.
