This paper presents a formal method for constructing solvable effusivity profiles, i.e. those leading to closed-form analytic one-dimensional solutions in the Laplace or Fourier domain for transient or steady-periodic temperature. A Liouville transformation applied to the heat equation expressed either in temperature or in heat flux yields two stationary Schrodinger-type equations. They reveal that temperature and heat flux depend solely on the thermal effusivity profile expressed in terms of the heat diffusion-time variable. Each Schrodinger-type equation can be cast in a system of two homologous Schrodinger-type equations, it means equations with the same potential function. The Darboux transformation is known to change a solvable linear second order differential equation into another solvable, but more complex, equation of the same type. We show that when applied to the heat transfer equation for continuously heterogeneous media and starting from an elementary solution, it allows finding simultaneously a set of more complex solvable effusivity profiles and the associated temperature solutions. The whole procedure constitutes a joint Property & Field Darboux Transformation (PROFIDT). By iterating these transformations we may construct the temperature solution related to an arbitrarily complex effusivity profile. Another way consists in applying an n-order Darboux-Crum transformation simultaneously on a seed field-function and on its associated square root effusivity profile (or the reciprocal). We illustrate these procedures, starting from the elementary solutions associated to a linear, a hyperbolic or a trigonometric effusivity profile (or their reciprocal). The solvable profiles and the temperature obtained therefrom present the distinct advantage of involving only elementary functions. As an alternative to the n-order solutions, we introduce the concept of solvable splines. They are obtained by stitching profiles from a set of parsimonious basis functions obtained after the application of a single PROFIDT. The resulting interpolation spline may be C-2 continuous at the nodes. The associated exact temperature solution is then obtained by the thermal quadrupole method (analytical transfer matrix method). The proposed methodology can be applied to a variety of physical problems where the dynamic field inside the graded material is described by a diffusion-like or a wave-like equation. (C) 2016 Elsevier Ltd. All rights reserved.