I'd like to give an alternative derivation not involving the clever (mystifying?) transformation to the heat equation and thus present a more general technique for solving constant coefficeint advection-diffusion PDEs. All we need is the Fourier transform:
\begin{align*}
\mathcal{F}[f] & = \int_{-\infty}^\infty e^{-i \omega y} f(y) dy,
\end{align*}
where $\mathcal{F}[f] = \mathcal{F}[f](\omega)$.
We'll use the following well-known facts:
\begin{align*}
\mathcal{F}\left[ \frac{1}{s \sqrt{2\pi}} \exp\left( -\frac{1}{2}\left(\frac{y - m}{s}\right)^2\right) \right] & = \exp(-i \omega m - s^2\omega^2/2), \\
\mathcal{F}\left[\frac{ \partial^n f}{\partial y^n}\right] & = (i \omega)^n \mathcal{F}[f], \\
\mathcal{F}[cf] & = c \mathcal{F}[f], \\
\mathcal{F}[f \ast g] & = \mathcal{F}[f]\mathcal{F}[g],
\end{align*}
where the convolution $(f \ast g)(y) = \int_{-\infty}^\infty f(z)g(y-z)dz$.
Using @olaker's notation, start with the constant coefficient PDE ((3) above)
$$
\frac{\partial C}{\partial \tau} = \frac{\sigma^2}{2}\frac{\partial^2 C}{\partial y^2} + \left(r-\frac{\sigma^2}{2}\right)\frac{\partial C}{\partial y} - rC.
$$
Step 1. Take the Fourier transform of each term term above and solve the resulting (very simple) ODE:
\begin{align*}
\frac{\partial \hat{C}}{\partial \tau} & = -\frac{\sigma^2 \omega^2}{2}\hat{C} + i \omega\left(r-\frac{\sigma^2}{2}\right)\hat{C} - r\hat{C}, \\
\hat{C} & = \hat{C}_0 e^{-r\tau} \exp \left(-\frac{\sigma^2 \omega^2}{2}\tau + i \omega\left(r-\frac{\sigma^2}{2}\right)\tau\right).
\end{align*}
Step 2. Letting $m = \left(\frac{\sigma^2}{2} - r\right)\tau$ and $s = \sigma \sqrt{\tau}$ from the Fourier transform notation, note
$$
\exp \left(-\frac{\sigma^2 \omega^2}{2}\tau + i \omega\left(r-\frac{\sigma^2}{2}\right)\tau\right) = \mathcal{F}\left[ \frac{1}{\sigma \sqrt{2\pi\tau}} \exp\left( -\frac{1}{2}\left(\frac{y - \left(\frac{\sigma^2}{2} - r\right)\tau}{\sigma \sqrt{\tau}}\right)^2\right) \right],
$$
so
\begin{align*}
\hat{C} & = \hat{C}_0 e^{-r\tau}\mathcal{F}\left[ \frac{1}{\sigma \sqrt{2\pi\tau}} \exp\left( -\frac{1}{2}\left(\frac{y - \left(\frac{\sigma^2}{2} - r\right)\tau}{\sigma \sqrt{\tau}}\right)^2\right) \right] \\
& = \frac{1}{\sigma \sqrt{2\pi\tau}} e^{-r\tau}\mathcal{F}\left[C_0 \ast \exp\left( -\frac{1}{2}\left(\frac{y - \left(\frac{\sigma^2}{2} - r\right)\tau}{\sigma \sqrt{\tau}}\right)^2\right) \right].
\end{align*}
Step 3. Take inverse transform:
\begin{align*}
C(y, \tau) & = \frac{1}{\sigma \sqrt{2\pi\tau}} e^{-r\tau} \int_{-\infty}^\infty C_0(z) \exp\left( -\frac{1}{2}\left(\frac{y - z - \left(\frac{\sigma^2}{2} - r\right)\tau}{\sigma \sqrt{\tau}}\right)^2\right) dz \\
C(y, \tau) & = \frac{1}{\sigma \sqrt{2\pi\tau}} e^{-r\tau} \int_{-\infty}^\infty C_0(z) \exp\left( -\frac{1}{2}\left(\frac{z - \left(y + \left(r - \frac{\sigma^2}{2}\right)\tau\right)}{\sigma \sqrt{\tau}}\right)^2\right) dz
\end{align*}
Step 4. Finally, change variables back $y \to S$, where we had $y = \log S$. Before we do this, note $S$ is really the ``initial'' stock price in the usual sense, i.e. $S = S_0$, but to be consistent we'll stick with $S$ as the initial (known) stock price. We'll also transform the $z$ variable, suggestively calling it $S_T$ by $S_T = e^z$.
\begin{align*}
C(S, \tau) & = \frac{1}{\sigma \sqrt{2\pi\tau}} e^{-r\tau} \int_0^\infty C_0(S_T) \frac{1}{S_T} \exp\left( -\frac{1}{2}\left(\frac{\log S_T - \left(\log S + \left(r - \frac{\sigma^2}{2}\right)\tau\right)}{\sigma \sqrt{\tau}}\right)^2\right) dS_T.
\end{align*}
Just note
\begin{align*}
f(S_T) := \frac{1}{S_T \sigma \sqrt{2\pi\tau}} \exp\left( -\frac{1}{2}\left(\frac{\log S_T - \left(\log S + \left(r - \frac{\sigma^2}{2}\right)\tau\right)}{\sigma \sqrt{\tau}}\right)^2\right)
\end{align*}
is the probability density function for a log $\mathcal{N}\left(\log S + \left(r - \frac{\sigma^2}{2}\right)\tau, \sigma^2 \tau\right)$ random variable, and under Black-Scholes, this is indeed the distribution of $S_T$ under $\mathbb{Q}$. Hence
\begin{align*}
C(S,\tau) = e^{-r\tau}\mathbb{E}_{\mathbb{Q}}[C_0(S_T)|\mathcal{F}_t].
\end{align*}
