  We present a novel numerical method for solving the anisotropic diffusion
equation in magnetic fields confined to a periodic box which is accurate and
provably stable. We derive energy estimates of the solution of the continuous
initial boundary value problem. A discrete formulation is presented using
operator splitting in time with the summation by parts finite difference
approximation of spatial derivatives for the perpendicular diffusion operator.
Weak penalty procedures are derived for implementing both boundary conditions
and parallel diffusion operator obtained by field line tracing. We prove that
the fully-discrete approximation is unconditionally stable. Discrete energy
estimates are shown to match the continuous energy estimate given the correct
choice of penalty parameters. A nonlinear penalty parameter is shown to provide
an effective method for tuning the parallel diffusion penalty and significantly
minimises rounding errors. Several numerical experiments, using manufactured
solutions, the ``NIMROD benchmark'' problem and a single island problem, are
presented to verify numerical accuracy, convergence, and asymptotic preserving
properties of the method. Finally, we present a magnetic field with chaotic
regions and islands and show the contours of the anisotropic diffusion equation
reproduce key features in the field.
