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 self-similar solutions of the associated self-similar PDE, we infer that the
t
−
1
/
2
t^{-1/2}
decay law of the diffusion amplitude is not necessary. In particular, it is possible to engineer setups of both the Cauchy problem and the initial-boundary value problem in which the solution decays at a different rate. Nevertheless, we observe that the
t
−
1
/
2
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 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.