I like to think of it in a similar way but with slightly different differential equations. (edit: below I managed to make it also intuitive for $\phi'(x) = -x \phi(x)$)
Case: heat equation
$$\frac{d}{dt} \left[ \frac{1}{\sqrt{t}} \phi(x/\sqrt{t}) \right]= 0.5 \frac{d^2}{dx^2} \left[ \frac{1}{\sqrt{t}} \phi(x/\sqrt{t}) \right]$$
This is how I actually like to see the normal distribution and the central limit theorem. The sum of variables which is in the limit like a diffusion process (which follows a differential equation)
For a random walk where the steps are some variable $X_i$ with unit variance and zero mean then you get that the distribution function is approaching a scaled normal distribution
$$P \left( \sum_{i=1}^t X_i = x\right) \approx \frac{1}{\sqrt{t}} \phi \left(x/\sqrt{t}\right) = f(x,t)$$
The change of this function can be seen as a bit simultaneously as the differential function for Brownian motion and is like a wave spreading out in time following the equation:
$$\frac{\partial}{\partial t} f(x,t) = \frac{\partial^2}{\partial x^2} f(x,t)$$
See also https://en.m.wikipedia.org/wiki/Normal_distribution#Exact_normality
Case: $\mathbf{\phi(x) +x\phi'(x)+ \phi''(x) =0}$
Now when we divide by $t$
$${P\left(\frac{\sum_{i=1}^t X_i }{t}= x\right) }\approx \phi\left({x}\right) = g(x,t)$$
and we could describe it as some scaled random walk like: if $Z_t = Z_{t-1} + X_t$ then for $Y_t = Z_t/\sqrt{t}$
$$Y_t -Y_{t-1} = - Y_{t-1} \frac{\sqrt{t-1}-\sqrt{t}}{\sqrt{t}} + X_t \frac{1}{\sqrt{t}}$$
which you could see as shrinking the current value by a factor and adding a variable $X_t$ scaled with some factor determined by $t$. Then in the limit the change of $\frac{\partial}{\partial t} g(x,t)$ should balance these two processes
$$ \frac{\partial g(x,t)}{\partial t} =
\underbrace{ a g(x,t) }_{\substack{
\text{shrinking}\\
\text{moves values up}
}}
+\underbrace{
\overbrace{x}^{
\substack{
\llap{\text{further }}
\rlap{\text{away the}} \\
\llap{\text{shrinking }}
\rlap{\text{is stronger}}\\ \,
}}
\frac{\partial g(x,t)}{\partial x}}_{
\substack{
\text{shrinking} \\
\text{shifts/squeezes the function}
}} +
\underbrace{ c \frac{\partial^2 g(x,t) }{\partial x^2 }}_{\text{diffusion}}=0$$
And the normal function is the function that such that the derivative of time (which is now expressed in terms of derivatives of space only) is zero.
So in this way we can relate intuitively $\phi(x) +x\phi'(x)+ \phi''(x) =0$ to a diffusion process with shrinking where the function remains unchanged.
Case: $\mathbf{x \phi(x) + \phi'(x) =0}$
This case is relatively similar to the one above.
The shrinking of the formula can be related to a flux that is pulling the density to the inside. The flux is the product of the amount of mass (which is moving) and the speed of the mass (which relates to the distance)
$$ \text{flux}_\text{shrinking} = -x \phi(x)$$
The diffusion can be related to a flux that is related to the slope of the function. If at some point there is more density to one direction than the other then the convolution/diffusion will cause some density to flow donwards on the slope.
$$ \text{flux}_\text{diffusion} = -\phi'(x)$$
When these two fluxes are opposite then there is no net flux and the function remains stable. So that is how you could view the relationship $\phi'(x) = -x\phi(x)$
I have made a computation to make the above intuitive idea more clear. In the computation I compute 1000 points according to some random distribution. And then I transform each point by scaling it with a factor $(1-c)$ and I am adding a centered Bernoulli variable to it with a factor $\sqrt{2c-c^2}$. This transformation will turn over time the distribution into a stable situation where the effect of the scaling is equal to the effect of the addition of the Bernoulli variable.
Below I have made two sketches for the intuition behind the terms in the differential equation.
(It is not a rigid derivation and one should go from the difference equations to differential and take the limit to very all linearizing the function, and also one could generalize the variable that is added, and represents the diffusion, which is now just a Bernoulli distributed variable. But I guess that in this way it is more intuitive and captures the essence more clearly)
# to plot points in the distribution
histpoints <- function(x, min, max) {
counts <- rep(0, length(min:max))
y <- rep(0,length(x))
for (i in 1:length(x)) {
bin <- ceiling(x[i]-min)
counts[bin] <- counts[bin]+1
y[i] <- counts[bin]
}
points(x,y, pch = 21, col = 1, bg = 1, cex = 0.4)
counts
}
# transforming the points by
# - scaling/shrinking
# - and adding a Bernoulli variable
convertpoints <- function(x,c,var = 25) {
x <- x * (1-c) # scaling
x <- x + sqrt(2*c-c^2) * (-1+2*rbinom(length(x), size = 1, prob = 0.5))*sqrt(var) # adding noise term
return(x)
}
# make 2000 points according to some funny distribution
set.seed(1)
start <- seq(-20,20,0.01)
x <- sample(start, 1000, replace = TRUE, prob = 20+start^2-(20^-2+20^-3)*start^4)
# plot initial histogram
layout(matrix(1:8,4))
par(mar=c(3,1,2,1))
hist(x, breaks = c(-40:40), xlim=c(-25,25), ylim = c(0,80), main = "begin", xlab = "", yaxt = "n", ylab = "", xaxt = "n")
bins <- histpoints(x,-30,30)
for (j in 1:7) {
for (i in 1:(100)) {
x <- convertpoints(x,0.003)
}
#plot transformed
hist(x, breaks = c(-40:40),
xlim=c(-25,25), ylim = c(0,80), main = paste0("after ",j*100," transforms"), xlab = "", yaxt = "n", ylab = "", xaxt = "n")
bins <- histpoints(x,-30,30)
}
