One starts with the Black-Scholes equation
$$\frac{\partial C}{\partial t}+\frac{1}{2}\sigma^2S^2\frac{\partial^2 C}{\partial S^2}+
rS\frac{\partial C}{\partial S}-rC=0,\qquad\qquad\qquad\qquad\qquad(1)$$
supplemented with the terminal and boundary conditions (in the case of a European call)
$$C(S,T)=\max(S-K,0),\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(2)$$
$$C(0,t)=0,\qquad C(S,t)\sim S\ \mbox{ as } S\to\infty.\qquad\qquad\qquad\qquad\qquad\qquad$$
The option value $C(S,t)$ is defined over the domain $0< S < \infty$, $0\leq t\leq T$.
Step 1. The equation can be rewritten in the equivalent form
$$\frac{\partial C}{\partial t}+\frac{1}{2}\sigma^2\left(S\frac{\partial }{\partial S}\right)^2C+\left(r-\frac{1}{2}\sigma^2\right)S\frac{\partial C}{\partial S}-rC=0.$$
The change of independent variables
$$S=e^y,\qquad t=T-\tau$$
results in
$$S\frac{\partial }{\partial S}\to\frac{\partial}{\partial y},\qquad \frac{\partial}{\partial t}\to - \frac{\partial}{\partial \tau},$$
so one gets the constant coefficient equation
$$\frac{\partial C}{\partial \tau}-\frac{1}{2}\sigma^2\frac{\partial^2 C}{\partial y^2}-\left(r-\frac{1}{2}\sigma^2\right)\frac{\partial C}{\partial y}+rC=0.\qquad\qquad\qquad(3)$$
Step 2. If we replace $C(y,\tau)$ in equation (3) with $u=e^{r\tau}C$, we will obtain that
$$\frac{\partial u}{\partial \tau}-\frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial y^2}-\left(r-\frac{1}{2}\sigma^2\right)\frac{\partial u}{\partial y}=0.$$
Step 3. Finally, the substitution $x=y+(r-\sigma^2/2)\tau$ allows us to eliminate the first order term and to reduce the preceding equation to the form
$$\frac{\partial u}{\partial \tau}=\frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial x^2}$$
which is the standard heat equation. The function $u(x,\tau)$ is defined for $-\infty < x < \infty$, $0\leq\tau\leq T$. The terminal condition (2) turns into the initial condition
$$u(x,0)=u_0(x)=\max(e^{\frac{1}{2}(a+1)x}-e^{\frac{1}{2}(a-1)x},0),$$
where $a=2r/\sigma^2$. The solution of the heat equation is given by the well known-formula
$$u(x,\tau)=\frac{1}{\sigma\sqrt{2\pi \tau}}\int_{-\infty}^{\infty} u_0(s)\exp\left(-\frac{(x-s)^2}{2\sigma^2 \tau}\right)ds.$$
Now if we evaluate the integral with our specific function $u_0$ and return to the old variables $(x,\tau,u)\to(S,t,C)$, we will arrive at the usual Black–Merton-Scholes formula
for the value of a European call. The details of the calculation can be found e.g. in The Mathematics of Financial Derivatives by Wilmott, Howison, and Dewynne (see Section 5.4).
