It is true that for any initial datum $u_0\in C^\infty(\mathbb{R})$ there exists a solution $u\in C^\infty(\mathbb{R}^+\times\mathbb{R})$ to the heat equation with initial condition $u(0,x)=u_0(x)$. As you point out, this will not be unique.
I can give a method of constructing such solutions now. The idea is to show that we can write $u=\sum_{n=1}^\infty f_n$ where $f_n$ are carefully constructed solutions chosen such that the partial sums $\sum_{n=0}^mf_n(0,x)$ eventually agree with $u_0(x)$ any bounded subset of the reals, and for which $f_n$ tends to zero arbitrarily quickly in the compact-open topology. First a bit of notation. I use $\mathbb{R}^+=[0,\infty)$ for the nonnegative reals. For a space $X$ then $C_0(X)$, $C^\infty(X)$, $C^\infty_0(X)$, and $C^\infty_K(X)$, represent the continuous real-valued functions on $X$ which are respectively vanishing at infinity, smooth, smooth and vanishing at infinity, and smooth with compact support. Let $(P_t)_{t\geq 0}$ be the kernels $$ \begin{align} &P_t\colon C_0(\mathbb{R})\to C_0(\mathbb{R}),\\\\ &P_tu(x)=\frac{1}{\sqrt{4\pi t}}\int_\mathbb{R}e^{-(x-y)^2/4t}u(y)\,dy \end{align} $$ for $t > 0$, and $P_0u=u$. This is the Markov transition function for Brownian motion (more precisely, for standard Brownian motion scaled by $\sqrt{2}$, because of the normalization used here). For $u\in C_K^\infty(\mathbb{R})$, then $f(t,x)=P_tu(x)$ is in $C_0^\infty(\mathbb{R}^+\times\mathbb{R})$, and is a solution to the heat equation with initial condition $f(0,x)=u(x)$, agreeing with the classical solution stated in the question. I'll also consider initial conditions $u\in C^\infty_K((a,\infty))$ (for $a\in\mathbb{R}$) by setting $u(x)\equiv0$ for all $x\le a$. The first step in the construction is to find initial conditions supported in $(a,\infty)$ so that $P_tu(0)$ approximates any given continuous function of time that we like.
- For any $a > 0$ and $h\in C_0((0,\infty])$, there exists a sequence $u_1,u_2,\ldots\in C^\infty_K((a,\infty))$ such that $\sqrt{t}P_tu_n(0)$ converges uniformly to $h(t)$ (over $t > 0$) as $n\to\infty$.
Consider the closure, $V$, in $C_0((0,\infty])$ (under the uniform norm) of the space of functions $t\mapsto\sqrt{4\pi t}P_tu(0)$ for $u\in C^\infty_K((a,\infty))$. Note that the limit as $t\to\infty$ does always exist and is just the integral of $u$. Then $V$ is a closed linear subspace of $C_0((0,\infty])$. Consider a sequence $u_n\in C_0((a,\infty))$ tending to the delta function $\delta_b$ at a point $b > a$, in the sense that $u_n$ all have support in the same compact set, and converge in distribution to $\delta_b$. Then, from the expression defining $P_t$, $\sqrt{4\pi t}P_tu_n(0)$ converges uniformly over $t > 0$ to $\exp(-b^2/4t)$. So, the function $t\mapsto\exp(-b^2/4t)$ is in $V$. As the set of functions of the form $t\mapsto\exp(-b^2/4t)$ for $b > a$ is closed under multiplication and separates points, the locally compact version of the Stone-Weierstrass theorem says that $V=C_0((0,\infty])$. Then, (1) follows.
- For any $u\in C^\infty_K((0,\infty))$ and $a,T > 0$, there exists a sequence $u_1,u_2,\ldots\in C^\infty_0((a,\infty))$ such that $f_n(t,x)\equiv P_t(u+u_n)(x)$ converges uniformly to zero (along with all its partial derivatives to all orders) over $t\in[0,T]$ and $x\le0$.
Choosing $0 < \epsilon < a$ so that the support of $u$ is contained in $(\epsilon,\infty)$, (1) implies that we can choose $u_n\in C^\infty_K((a,\infty)$ so that $\sqrt{t}P_tu_n(\epsilon)$ converges uniformly to $-\sqrt{t}P_tu(\epsilon)$ over $t\ge0$ as $n\to\infty$. Then, $f_n(t,x)\equiv P_t(u+u_n)(x)$ is a bounded solution to the heat equation with boundary conditions $f_n(0,x)=0$ for $x\le\epsilon$ and $f_n(t,\epsilon)=P_t(u+u_n)(\epsilon)$. It is then standard that the solution is given by an integral over the boundary, $$ f_n(t,x)=\int_0^t\frac{\epsilon-x}{\sqrt{4\pi (t-s)^3}}e^{-(\epsilon-x)^2/4(t-s)}\sqrt{s}P_s(u_n+u)(\epsilon)\frac{ds}{\sqrt{s}} $$ for $x\le0$. Differentiating this wrt $x$ and $t$ an arbitrary number of times, and using dominated convergence as $n\to\infty$, it follows that $f_n(t,x)$ (and all its partial derivatives) converge uniformly to zero over $x\le0$ and $t\in[0,T]$ as $n\to\infty$.
- Suppose that $T > 0$, $0 < a < b$ and $u\in C^\infty_K(\mathbb{R})$ has support contained in $(-\infty,-a)\cup(a,\infty)$. Then, there exists a sequence $u_n\in C^\infty_K(\mathbb{R})$ with supports in $(-\infty,-b)\cup(b,\infty)$ such that $P_t(u+u_n)(x)$ (and all its partial derivatives) tends to 0 uniformly over $\vert x\vert\le a$ and $t\in[0,T]$.
Set $u^+(x)=1_{\{x > 0\}}u(x)$ and $u^-(x)=1_{\{x < 0\}}u(x)$. Applying (2) to $u^+(x+a)$, there exists a sequence $u^+_n\in C^\infty_K(\mathbb{R})$ with supports in $(b,\infty)$ such that $P_t(u^++u^+_n)(x)$ tends to 0 uniformly over $x\le a$ and $t\in[0,T]$. Applying the same argument to $u^-(-x-a)$, there exists $u^-_n\in C^\infty_K(\mathbb{R})$ with support in $(-\infty,-b)$ such that $P_t(u^-+u^-_n)(x)$ tends to zero uniformly over $x\ge -a$ and $t\in[0,T]$. A sequence satisfying the statement of (3) is $u_n=u^+_n+u^-_n$.
- If $u\in C^\infty(\mathbb{R})$ then there exists $f\in C^\infty(\mathbb{R}^+\times\mathbb{R})$ solving the heat equation with initial condition $f(0,x)=u(x)$.
We can inductively choose a sequence $u_n\in C^\infty_K(\mathbb{R})$ such that $\sum_{m=1}^nu_m(x)=u(x)$ for $\vert x\vert\le n+1$ and $n\ge1$. Let $u_1=u$ on $[-2,2]$ and then, for each $n\ge2$, apply the following step.
- As $\tilde u=u-\sum_{m=1}^{n-1}u_m$ is zero on $[-n,n]$, it has support in $(-\infty,-a)\cup(a,\infty)$ for $a=n-1/2$. Choosing $b=n+1$, by (3), there exists $v\in C^\infty_K(\mathbb{R})$ with support in $(-\infty,-b)\cup(b,\infty)$ such that $P_t(\tilde u +v)(x)$ along with all its partial derivatives up to order $n$ are bounded by $2^{-n}$ over $\vert x\vert\le n-1/2$ and $t\le n$. Take $u_n=\tilde u + v$.
Setting $f_n(t,x)=P_tu_n(x)$, then $f_n$ are smooth functions satisfying the heat equation and the initial conditions $\sum_{m=1}^nf_m(0,x)=u(x)$ for $\vert x\vert\le n+1$. Also, by the choice of $u_n$, for $n > 1$ then $f_n$ together with all its derivatives up to order $n$ is bounded by $2^{-n}$ on $[0,n]\times[1-n,n-1]$. As $\sum_n2^{-n} < \infty$, the limit $f=\sum_nf_n$ exists and is smooth with all partial derivatives commuting with the summation. Then $f$ has the required properties.
