  The diffusion equation is a universal and standard textbook model for partial
differential equations (PDEs). In this work, we revisit its solutions, seeking,
in particular, self-similar profiles. This problem connects to the classical
theory of special functions and, more specifically, to the Hermite as well as
the Kummer hypergeometric functions. Reconstructing the solution of the
original diffusion model from novel self-similar solutions of the associated
self-similar PDE, we infer that the $t^{-1/2}$ decay law of the diffusion
amplitude is {\it not necessary}. In particular, it is possible to engineer
setups of {\it both} the Cauchy problem and the initial-boundary value problem
in which the solution decays at a {\it different rate}. Nevertheless, we
observe that the $t^{-1/2}$ rate corresponds to the dominant decay mode among
integrable initial data, i.e., ones corresponding to finite mass. Hence, unless
the projection to such a mode is eliminated, generically this decay will be the
slowest one observed. In initial-boundary value problems, an additional issue
that arises is whether the boundary data are \textit{consonant} with the
initial data; namely, whether the boundary data agree at all times with the
solution of the Cauchy problem associated with the same initial data, when this
solution is evaluated at the boundary of the domain. In that case, the power
law dictated by the solution of the Cauchy problem will be selected. On the
other hand, in the non-consonant cases a decomposition of the problem into a
self-similar and a non-self-similar one is seen to be beneficial in obtaining a
systematic understanding of the resulting solution.
