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 χ τ ∼ k −α where χ is the parallel diffusivity, and α = 1 (α = 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 = (χ t) −γ/2 L (χ t) −γ/2 δψ , where δψ is a radial coordinate. In the local case, f is Gaussian and the scaling is sub-diffusive, γ = 1/2. In the non-local case, f decays algebraically, L(η) ∼ η −3 , and the scaling is diffusive, γ = 1.