(This is a follow-up to Gjergji Zaimi's beautiful answer: a minor simplification of the method applied there rather than a different solution. Written as a separate answer upon request of the original poster.)

Let $g_s(x) = (\pi s)^{-1/2} \exp(x^2 / s)$ denote the Gauss–Weierstrass kernel. Note that $g_s * g_t = g_{s + t}$.

Choose $\delta \in (0, 1)$ sufficiently small, so that
$$ \begin{aligned} h(x) & = \sum_{i = 1}^n a_i \delta^{-1/2} \exp((x - x_i)^2 / \delta) - c \\ & = \pi^{1/2} \sum_{i = 1}^n a_i g_\delta(x - x_i) - c \end{aligned} $$
has exactly $2 n$ zeroes; or, more precisely, it changes its sign exactly $2 n$ times. (Is is intuitively clear, and relatively straightforward to prove, that $\delta > 0$ small enough have this property, but the details are a bit messy.) Then $f$ can be written as a convolution of $h$ and $g_{1 - \delta}$:
$$ \begin{aligned} f(x) & = \sum_{i = 1}^n a_i \exp((x - x_i)^2) - c \\ & = \pi^{1/2} \sum_{i = 1}^n a_i g_1(x - x_i) - c \\ & = \biggl(\pi^{1/2} \sum_{i = 1}^n a_i g_\delta(\cdot - x_i) - c\biggr) * g_{1 - \delta}(x) . \end{aligned} $$
However, convolution with the Gauss–Weierstrass kernel does not increase the number of sign-changes: $g_s$ is a *Pólya frequency function*, and hence a *variation diminishing kernel*. For more on this, see, for example, review sections written by Samuel Karlin in *Selected works of I. J. Schoenberg*. It follows that $f$ changes its sign at most $2 n$ times, as claimed.