  A novel approach that enables the study of parallel transport in magnetized
plasmas is presented. The method applies to general magnetic fields with local
or nonlocal parallel closures. Temperature flattening in magnetic islands is
accurately computed. For a wave number $k$, the fattening time scales as
$\chi_{\parallel} \tau \sim k^{-\alpha}$ where $\chi$ is the parallel
diffusivity, and $\alpha=1$ ($\alpha=2$) for non-local (local) transport. The
fractal structure of the devil staircase temperature radial profile in weakly
chaotic fields is resolved. In fully chaotic fields, the temperature exhibits
self-similar evolution of the form $T=(\chi_{\parallel} t)^{-\gamma/2} L \left[
(\chi_{\parallel} t)^{-\gamma/2} \delta \psi \right]$, where $\delta \psi$ is a
radial coordinate. In the local case, $f$ is Gaussian and the scaling is
sub-diffusive, $\gamma=1/2$. In the non-local case, $f$ decays algebraically,
$L (\eta) \sim \eta^{-3}$, and the scaling is diffusive, $\gamma=1$.
