In magnetic toroidal confinement fusion, transport of particles or temperature happens primarily along magnetic field lines, with the ratio of diffusion parallel and perpendicular to the field lines sometimes exceeding the order of $10^8$. A simple model for the transport is the anisotropic diffusion equation, however it can be difficult to solve numerically due to the ratio of diffusion coefficients, since the numerical error can quickly overwhelm the diffusion perpendicular to field lines. In this presentation we discuss a novel penalty approach for computing the solution to the anisotropic diffusion equation in periodic geometry which is provably stable. We use summation by parts finite differences to compute the slow scale diffusion. Fast scale diffusion is computed using interpolation and applied with a non-linear penalty term which is determined for numerical stability. Recent improvements to our method have extended our method from Cartesian to cylindrical geometry, allowing us to compute solutions in geometry more closely resembling those involved in fusion plasmas. We present new results in slab geometry showing how the grid can designed to efficiently resolve regions of interest. Further we show results in cylindrical geometry and show how the method can be extended to the geometry of confinement fusion devices.