  This paper presents a detailed discussion of the ``Newton's method''
algorithm for finding apparent horizons in 3+1 numerical relativity.
  We describe a method for computing the Jacobian matrix of the finite
differenced $H(h)$ function by symbolically differentiating the finite
difference equations, giving the Jacobian elements directly in terms of the
finite difference molecule coefficients used in computing $H(h)$. Assuming the
finite differencing scheme commutes with linearization, we show how the
Jacobian elements may be computed by first linearizing the continuum $H(h)$
equations, then finite differencing the linearized (continuum) equations.
  We find this symbolic differentiation method of computing the $H(h)$ Jacobian
to be {\em much} more efficient than the usual numerical perturbation method,
and also much easier to implement than is commonly thought.
  When solving the discrete $H(h) = 0$ equations, we find that Newton's method
generally converges very rapidly. However, if the initial guess for the horizon
position contains significant high-spatial-frequency error components, Newton's
method has a small (poor) radius of convergence. This is {\em not} an artifact
of insufficient resolution in the finite difference grid; rather, it appears to
be caused by a strong nonlinearity in the continuum $H(h)$ function for
high-spatial-frequency error components in $h$. Robust variants of Newton's
method can boost the radius of convergence by O(1) factors, but the underlying
nonlinearity remains, and appears to worsen rapidly with increasing
initial-guess-error spatial frequency.
  Using 4th~order finite differencing, we find typical accuracies for computed
horizon positions in the $10^{-5}$ range for $\Delta\theta = \frac{\pi/2}{50}$.
