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.
1) 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.
2) 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$.
3) 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$.
4) 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.
Best Answer
(Not sure if I understand the question correctly.)
If $p_{t,x}(y) = p(t, x, y)$ is the fundamental solution (a.k.a. the heat kernel), then $p_{t,x}$ converges as $t \to 0^+$ to the Dirac measure $\delta_0$ in the sense of weak* convergence of measures, and hence also in the sense of distributions.
For any initial value given by a distribution $f$ in $(-1,1)$ (say: compactly supported, but this can be extended slightly), $u(t, x) = \langle f, p_{t,x}$ makes sense. Then $u$ can be proved to solve the heat equation, and $u(t, \cdot)$ converges in the space of distributions to $f$ as $t \to 0^+$.
I do not know the literature well, but I would check any book on applications of the distribution theory to PDEs for more on that. For example, Section 16.7 in Vladimirov's Methods of the theory of generalized functions might be a good reference.