I'm trying to evaluate the integral $$ \Psi(\mu,\nu)=\int_{-\infty}^{\infty}\operatorname{e}^{-\mu x^2}f(\nu x)\operatorname{d}\!x\qquad(\text{for}\; \mu>0)\tag 1 $$ where $\nu\in\Bbb R$ and $f$ satisfies the functional equation $$ f(x+y)+f(x-y)=2f(x)f(y).\tag 2 $$
I tried to evaluate the product $\Psi(\mu,\nu)\cdot \Psi(\mu,\rho)$ and, after some manipulations, I found $$ \Psi(\mu,\nu)\cdot \Psi(\mu,\rho)=\gamma \Psi\left(\mu,\sqrt{\nu^2+\rho^2}\right)\tag 3 $$ where $\gamma=\sqrt{\pi/\mu}$. Putting $\Phi(\alpha)=\frac{1}{\gamma}\Psi(\mu,\alpha)$, the equation (3) becomes $$ \Phi(\nu)\Phi(\rho)=\Phi\left(\sqrt{\nu^2+\rho^2}\right).\tag 4 $$
Do you have any suggestion on how to procede?
To avoid confusion, I'm using $\Phi_\mu$ to represent $\Phi$ throughout. If $f(0) \neq 0$ and $\Phi_\mu$ is integrable on finite intervals for every $\mu$, then $\Phi_\mu(\nu)$ is of the form $e^{c(\mu)\nu^2}$ for some function $c(\mu)$. To see this, fix $\mu > 0$. Set $\nu = \rho = 0$ in the functional equation for $\Phi_\mu$ to get $\Phi_\mu(0)^2 = \Phi_\mu(0)$. Since $f(0) \neq 0$, $\Phi_\mu(0) \neq 0$, so $\Phi_\mu(0) = 1$. Hence, setting $\rho = 0$ in the functional equation for $\Phi_\mu$ results in $\Phi_\mu(\nu) = \Phi_\mu(|\nu|)$ for all $\nu \in \Bbb R$. Consider the function $g_\mu$ by the equation $g_\mu(\nu) = \Phi_\mu(\sqrt{\nu})$, for all $\nu \ge 0$. Then $g_\mu$ satisfies the functional equation $$g_\mu(\nu) g_\mu(\rho) = g_\mu(\nu + \rho).$$ Since $\Phi_\mu$ is integrable on finite intervals, so is $g_\mu$. Let $\nu_0$ be a point such that $\int_0^{\nu_0} g_\mu(\rho)\, d\rho \neq 0$. Then $$g_\mu(\nu) = \frac{\int_0^{\nu_0} g_\mu(\nu)g_\mu(\rho)\, d\rho}{\int_0^{\nu_0} g_\mu(\rho)\, d\rho} = \frac{\int_0^{\nu_0} g_\mu(\nu + \rho)\, d\rho}{\int_0^{\nu_0} g_\mu(\rho)\, d\rho} = \frac{\int_{\nu}^{\nu + \nu_0} g_\mu(\rho)\, d\rho}{\int_0^{\nu_0} g_\mu(\rho)\, d\rho}.$$ So $g_\mu$ is continuous and hence differentiable. Let $c(\mu) := g_\mu'(0)$. Then \begin{align}g_\mu'(\nu) &= \lim_{h\to 0} \frac{g_\mu(\nu + h) - g_\mu(\nu)}{h}\\ &= \lim_{h\to 0} \frac{g_\mu(\nu)g_\mu(h) - g_\mu(\nu)}{h}\\ & = g_\mu(\nu) \lim_{h\to 0} \frac{g_\mu(h) - 1}{h}\\ &= g_\mu(\nu) \lim_{h\to 0} \frac{g_\nu(h) - g_\nu(0)}{h}\\ & = c(\mu)g_\mu(\nu).\end{align} Therefore $g_\mu(\nu) = e^{c(\mu)\nu}$, and $\Phi_\mu(\nu) = g_\mu(\nu^2) = e^{c(\mu)\nu^2}$ for all $\nu \ge 0$. Since $\Phi_\mu(\nu) = \Phi_\mu(|\nu|)$ for all $\nu \in \Bbb R$, $\Phi_\mu(\nu) = e^{c(\mu)\nu^2}$ for all $\nu \in \Bbb R$.