$\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\om}{\omega}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}}
\renewcommand{\P}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}
\newcommand{\tf}{\widetilde{f}}$
A brief answer, at least to the "Math (integro-differential equations)" question, is no. For example, let $U(x)=\cosh ax$ for $a\ne0$. Then
\begin{equation}
\frac{dEU(x+N(0,1))}{dx}=c\ \frac{dU(x)}{dx} \tag{1}
\end{equation}
for $c=e^{a^2/2}$, but
$U' U''' \ne U'' U''$.
A (much) longer story is to consider a general solution of equation (1), which can be rewritten as
\begin{equation*}
V*\phi=cV,
\end{equation*}
where $V:=U'$ and $\phi$ is the standard normal density. Formally, considering the Fourier transform $\hat V(t):=\int_\R e^{itx}V(x)\,dx$, we have
\begin{equation}
(\sqrt{2\pi}\,\phi-c)\hat V=0 \tag{1.5}
\end{equation}
and "hence"
\begin{equation}
\hat V=\sum_j C_j\de_{iu_j}, \tag{2}
\end{equation}
where the $C_j$'s are constants, the $iu_j$'s are the complex roots of the equation $\sqrt{2\pi}\,\phi-c=0$, and $\de_b$ is the Dirac distribution at point $b$, so that $\langle \de_b,f\rangle=f(b)$ for test functions $f$. "Therefore",
$V(x)$ is a linear combination of the $e^{u_jx}$'s and hence $U(x)$ is a linear combination of the $e^{u_jx}$'s and of a constant (and also of $x$ if $c=1$). This is how the solution $U(x)=\cosh ax$ was found.
On the other hand, it is easy to verify by direct calculation that any linear combination of the $e^{u_jx}$'s (with the $u_j$'s as above) is a solution to equation (1). Particular examples of such linear combinations are $U(x)=\cosh ax$ and $U(x)=\sinh ax$.
Remaining problems with the above approach to finding a general solution of equation (1) are as follows:
We need to define an appropriate class of functions $U$ for which equation (1) is considered.
Since $\hat V$ will in general be a distribution, we need to define an appropriate space of test functions. The space of all compactly supported smooth functions will not do, since the Fourier transform does not map this space into itself. The space of all tempered test functions will not do here either, since we should allow $U$ (and hence $V$) to grow at least exponentially fast. Perhaps the space of all functions $f$ such that
\begin{equation}
|f(x)|\le Ce^{-ax^2}\quad\text{and}\quad |\hat f(x)|\le Ce^{-bx^2} \tag{3}
\end{equation}
for some positive $C,a,b$ (possibly depending on $f$) and all real $x$ will do as the test functions, but
(i) this will exclude functions such as $V(x)=e^{x^2/4}$, for which $V*\phi$ is well defined and
(ii) I am not aware of a Fourier transform theory for the distributions corresponding to the test functions as in (3); cf. Hardy, A Theorem Concerning Fourier Transforms. Here we also to need to rigorously justify the transition from (1.5) to (2).
- We need to attach a rigorous meaning to the linear combination in (2).
