What you call "combination of variables" I think is often called "separation of variables" in my humble sphere of math. Let's take an example:
$$u_{xx}+u=u_t$$
Here, $u$ is the dependent variable and $x,t$ are independent. Usually, we suppose $u(x,t) = X(x)T(t)$ hence $u_{xx} = X''T$ and $u_t = XT'$. Substituting into the PDE gives,
$$X''T+XT=XT' $$
If the solution we seek is nontrivial then division by $XT$ is allowed,
$$ \frac{X''}{X}+1 = \frac{T'}{T} $$
At this point we note that the left and right sides of this equation are functions
of variables which we took as independent. It follows there is some constant $K$ such that:
$$ \frac{X''}{X}+1 = \frac{T'}{T} =K$$
Hence, we find solving the PDE reduces to solving two ODEs:
$$ T' = KT \qquad X''+(1-K)X=0. $$
Notice that $K$ is arbitrary with the data given. This means that we don't
really have just two ODES, we have infinitely many. What typically happens is that
in a well-posed problem you are also given Boundary Values (BV) and an initial condition.
In nice problems the BVs will force $K$ to be a certain discrete set of values, countably many in fact. Then for each choice of $K$ you get two ODEs and a couple of solutions to the common BV. After all of this you sum these together to form the formal solution and match it to the initial condition(s) by the Fourier technique. So, in that sense there are not just a couple constants as in the solution to a single ODE. There is infinite flexibility in the solution because you have a family of ODEs forming the solution.
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
$$x^2\dfrac{\partial u}{\partial x}+y\dfrac{\partial u}{\partial y}=1-xyu$$ The Charpit-Lagrange characteristic ODEs are : $$\frac{dx}{x^2}=\frac{dy}{y}=\frac{du}{1-xyu}$$ A first characteristic equation comes from solving $\frac{dx}{x^2}=\frac{dy}{y}$ : $$y\:e^{1/x}=c_1$$ A second characteristic equation comes from solving $\frac{dx}{x^2}=\frac{du}{1-xyu}\quad\implies\quad \frac{dx}{x^2}=\frac{du}{1-x(c_1e^{-1/x})u}$
$\frac{du}{dx}+c_1\frac{1}{x}e^{-1/x}u=\frac{1}{x^2}\quad$ is a first order linear ODE which solution is : $$u\:e^{-c_1\text{Ei}(-1/x)}-\int_1^x \left(e^{c_1\text{Ei}(-1/\xi)}\right)\frac{d\xi}{\xi^2}=c_2$$ Ei is the exponential integral.
The general solution of the PDE on the form of implicit equation $c_2=F(c_1)$ is : $$u\:e^{-c_1\text{Ei}(-1/x)}-\int_1^x \left(e^{c_1\text{Ei}(-1/\xi)}\right)\frac{d\xi}{\xi^2}=F(c_1)$$ $$u=e^{c_1\text{Ei}(-1/x)}\left(\int_1^x \left(e^{c_1\text{Ei}(-1/\xi)}\right)\frac{d\xi}{\xi^2}+F(c_1)\right)$$ F is an arbitrary function, to be determined according to some boundary condition (missing in the wording of the question). $$\boxed{u(x,y)=e^{y\:e^{1/x}\text{Ei}(-1/x)}\left(\int_1^x \left(e^{y\:e^{1/x}\text{Ei}(-1/\xi)}\right)\frac{d\xi}{\xi^2}+F(y\:e^{1/x})\right)}$$