For your question 1, I prefer to use the substitution $\eta=\frac{y}{\sqrt{\alpha t}}$ first, then substitute $f$ into the PDE.
$$u\Big(\frac{y}{\sqrt{\alpha t}},1\Big) = f\Big(\frac{y}{\sqrt{\alpha t}}\Big)$$
Say $\eta=\frac{y}{\sqrt{\alpha t}}$, so $u=f(\eta)$.
For use later:
\begin{align}
\frac{\partial \eta}{\partial t}=&-\frac{1}{2}\frac{y}{\sqrt{\alpha}}t^{-3/2}\\
\frac{\partial \eta}{\partial y}=&\frac{1}{\sqrt{\alpha t}}\\
\frac{\partial y^*}{\partial y}=&\frac{1}{\sqrt{\alpha t_0}}\\
\frac{\partial t^*}{\partial t}=&\frac{1}{\sqrt{t_0}}.
\end{align}
For the LHS of the PDE, using the chain rule (and also using the condition $t_0=t$):
\begin{align}
\frac{\partial u}{\partial t^*}=&\frac{\partial f(\eta)}{\partial \eta}\frac{\partial \eta}{\partial t}\frac{\partial t}{\partial t^*}\\
=&f'(\eta)\Big(-\frac{1}{2}\frac{y}{\sqrt{\alpha}}t^{-3/2}\Big)t_0\\
=&f'(\eta)\Big(-\frac{1}{2}\frac{y}{\sqrt{\alpha}}t^{-3/2}\Big)t\\
=&-\frac{1}{2}f'(\eta)\Big(\frac{y}{\sqrt{\alpha t}}\Big)\\
=&-\frac{1}{2}\eta f'(\eta).
\end{align}
For the RHS of the PDE, a similar process is used:
\begin{align}
\frac{\partial}{\partial y^*}\Big(\frac{\partial u}{\partial y^*}\Big)=&\frac{\partial}{\partial \eta}\Big(\frac{\partial f(\eta)}{\partial \eta}\frac{\partial \eta}{\partial y}\frac{\partial y}{\partial y^*}\Big)\frac{\partial \eta}{\partial y}\frac{\partial y}{\partial y^*}\\
=&\frac{\partial}{\partial \eta}\Big(f'(\eta)\frac{1}{\sqrt{\alpha t}}\sqrt{\alpha t_0}\Big)\frac{1}{\sqrt{\alpha t}}\sqrt{\alpha t_0}\\
=&f''(\eta).
\end{align}
Therefore
$$-\frac{1}{2}\eta f'(\eta)=f''(\eta)$$
which if you want, you can substitute $\eta$ back in to produce
$$-\frac{1}{2}\frac{y}{\sqrt{\alpha t}} f'\Big(\frac{y}{\sqrt{\alpha t}}\Big)=f''\Big(\frac{y}{\sqrt{\alpha t}}\Big)$$
as you required.
Continuing on for your question 2 for completeness, the other commenters and answer are right. Use
$$\Theta(\eta)=\frac{\partial f(\eta)}{\partial \eta}$$
to change the PDE to
$$\frac{\partial \Theta(\eta)}{\partial\eta}=-\frac{1}{2}\eta \Theta(\eta)$$
which has the solution of
$$\Theta(\eta)=Ce^{-\eta^2/4}.$$
Integrate this
$$f(\eta)=C\int_\eta^\infty e^{-\eta^2/4}d\eta$$
then change variables using $\xi^2=\eta^2/4$ and $\eta=2d\xi$ to get
$$f(\xi)=2C\int_{\eta/2}^\infty e^{-\xi^2}d\xi. \tag{$\star$}$$
By definition (good example here), the error function of parameter $z$ is
$$\text{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^z e^{-\phi^2}d\phi$$
and has the properties $\text{erf}(-z)=-\text{erf}(z)$, $\text{erf}(0)=0$, and $\text{erf}(\infty)=1$.
Using these properties, equation ($\star$) can be split up into
$$f(\xi)=2C\Bigg[\int_{0}^\infty e^{-\xi^2}d\xi-\int_{0}^{\eta/2} e^{-\xi^2}d\xi\Bigg].$$
The first integral evaluates to $\sqrt{\pi}/2\text{erf}(\infty)$, and the second integral evaluates to $\sqrt{\pi}/2\text{erf}(\eta/2)$. Using the definition of the error function complement
\begin{align}
f(\xi)=&2C\frac{\sqrt{\pi}}{2}\Big[1-\text{erf}(\eta/2)\Big]\\
=&C\sqrt{\pi}\text{erfc}(\eta/2)
\end{align}
Then apply your boundary conditions and find constant $C$. I attempted this but got some unexpected results. I believe your boundary conditions may be incorrect; if $f\rightarrow 1$ as $\eta\rightarrow \infty$, that means the value $u$ that is diffusing is non-zero at $y=\infty$. A much more reasonable condition is $f\rightarrow 0$ as $\eta\rightarrow \infty$. In addition the other boundary condition may be incorrect; $f(0)=0$ means that $u=0$ when $\eta=0$, meaning at the initial time and origin, there is no value for $u$. I'm not sure this is realistic for the system you're attempting to model. Perhaps check the boundary conditions and come back.
Finding substitution using dimensional analysis
As I mentioned in a comment, It's actually quite interesting how the substitution $\eta=\frac{y}{\sqrt{\alpha t}}$ is chosen. Most of the time is seems they are just pulled out of someone's imagination, but you can derive it using dimensional analysis.
This is part of a method called the Buckingham Pi theorem, which is explained really well in the book "Scaling" by G.I Barenblatt. The first part of the work is already done, where you choose a non-dimensionalisation factor for $T$. I'll call it $\Pi$ where you called it $u$, to be consistent with the original theorem
$$ \Pi = \frac{T-T_w}{T_i-T_w}. $$
Then list the other variables in the system and note down their dimensions. For this I'm assuming $T$ is some scalar, like temperature (with dimension symbol $\Theta$).
\begin{align}
[T]=&\Theta\\
[T_w]=&\Theta\\
[T_i]=&\Theta\\
[\alpha]=&L^2T^{-1}\\
[y]=&L\\
[t]=&T
\end{align}
Note that (when grouping all the temperatures together) there are 4 variables but only 3 units between them: $\Theta,T,L$. Now partition the variables into two groups, variables with independent dimensions, and variables who's dimensions can be expressed in terms of the dimensions of the variables with independent dimensions. In this case, the variables with independent dimensions can be chosen to be $T,\alpha,t$.
If we look at $y$, which is one of our variables who's dimensions can be expressed in terms of the dimensions of the variables with independent dimensions, we can see that
$$[y]=[T]^{p_0}[\alpha]^{p_1}[t]^{p_2}.$$
Through deduction, $p_0=0$, $p_1=\frac{1}{2}$, and $p_2=\frac{1}{2}$. Therefore
$$[y]=[\alpha]^{1/2}[t]^{1/2}$$
and we can define a non-dimensional parameter
$$\Pi_1=\frac{y}{\alpha^{1/2}t^{1/2}}.$$
Then using the Buckingham Pi theorem (quoted from "Scaling", Barenblatt, 2003):
... a physical relationship between some dimensional (generally speaking) quantity and several dimensional governing parameters can be rewritten as a relationship between a dimensionless parameter and several dimensionless products of the governing parameters; the number of dimensionless products is equal to the total number of governing parameters minus the number of governing parameters with independent dimensions.
Recall that we have 4 variables, and 3 variables with independent dimensions. But since we use the non-dimensional version of $T$ this removes it from the relationship (from $4$ to $3$), so $3-3=0$ dimensionless products are needed to represent the relationship. So $\Pi$ is written as purely the dimensionless parameter, which is a function of $\Pi_1$
$$\Pi=f(\Pi_1).$$
Going back to our original notation, the substitution parameter is revealed
$$u=f\Big(\frac{y}{\sqrt{\alpha t}}\Big).$$
I hope that was understandable. The actual proof of the Buckingham Pi theorem is not the hardest to understand, but it's pretty long. I really recommend reading "Scaling" for the proof, especially if you're studying physics. It's amazing how it's possible to reduce the complexity of physical relations to only a few parameters.
Best Answer
It is perfectly possible, to show you why use the chain rule. Here's (perhaps) a simpler version, call
\begin{eqnarray} T &=& (T_i - T_w)u \\ r &=& a r^* \\ t &=& t_0 t^* \end{eqnarray}
So that
$$ \frac{\partial }{\partial r} = \frac{{\rm d}r^*}{{\rm d}r}\frac{\partial }{\partial r^*} = \frac{{\rm d}(r/a)}{{\rm d}r}\frac{\partial}{\partial r^*} = \frac{1}{a}\frac{\partial}{\partial r^*} $$
Similarly
$$ \frac{\partial}{\partial t} = \frac{1}{t_0}\frac{\partial}{\partial t^*} $$
And your equation becomes
\begin{eqnarray} \require{cancel} \frac{1}{t_0}\frac{\partial}{\partial t^*}[\color{blue}{\cancel{(T_i - T_w)}}u] &=& \frac{\alpha}{a^2}\frac{1}{r^*}\frac{\partial }{\partial r^*}\left( \frac{\color{red}{\cancel{a}} r^*}{\color{red}{\cancel{a}}} \frac{\partial }{\partial r^*}[\color{blue}{\cancel{(T_i - T_w)}}u]\right) \\ \implies~~\frac{1}{t_0} \frac{\partial u}{\partial t^*} &=& \frac{\alpha}{a^2} \frac{1}{r^*}\frac{\partial}{\partial r^*}\left(r^* \frac{\partial u}{\partial r^*} \right) \end{eqnarray}