Finding the solution of $u_x + y u_y = 0$ using $u(0, y) = y^3$

You seem to make some confusion due to the use of the variables $x,y$ for multiple (and different) purposes. The auxiliary (initial value) equation you write as $u(0,y)=y^3$ but you could equally well write it as $u(0,t)=t^3$ for any real value of $t$ (better not use $y$ here).

For a given point $(x,y)$ in the plane you have from the result just before (7): $u(x,y) = u(0,ye^{-x})$. Setting $t=y e^{-x}$ and using the previous we get:

$$ u(x,y)=u(0,ye^{-x}) = u(0,t) = t^3 = (y e^{-x})^3 = y^3 e^{-3x}.$$

Hope this brings you some peace of mind.


Let us verify the solution, using the Cauchy method. Trying to find the solution in the form of $$u(x,y) = X(x)Y(y),\tag1$$

one can get $$X'(x)Y(y) + yX(x)Y'(y) = 0,$$ $$\frac{X'(x)}{X(x)} + y\frac{Y'(y)}{Y(y)} = 0.\tag2$$

Since the first term of $LHS(2)$ does not depend of $y$ and the second term does not depend of $x,$ then solutions of $(2)$ exists only if these terms are the opposite constants.

Assume the second term as $\lambda,$ then $$\frac{X'(x)}{X(x)} = -\lambda,\quad \frac{Y'(y)}{Y(y)} = \frac\lambda y,\tag3$$ $$ \begin{cases} \ln C_1^{-1}X = - \lambda x\\[4pt] \ln C_2^{-1}Y = \lambda \ln y \end{cases}\Rightarrow \begin{cases} X = C_1e^{-\lambda x}\\[4pt] Y = C_2y^\lambda, \end{cases} $$

$$u(x,y,\lambda) = C_\lambda(y e^{-x})^\lambda.\tag4$$

The general solution can be defined as the arbitrary linear combination of such solutions over the domain $D(\lambda)$ in the form of $$u(x,y) = \int\limits_{D(\lambda)} w(\lambda)\ y^\lambda\ e^{-\lambda x}\text{ d}\lambda = f(e^{-x}y).\tag5$$

If $u(0,y) = y^3,$ then

  • from $(4)$ should $\lambda = 3, C_\lambda = 1,u(x,y,3) = y^3e^{-3x};$
  • from $(5)$ should $w(\lambda) = \delta(\lambda-3),\quad u(x,y)= y^3e^{-3x}.$

Also, some details by the paper.

  • The directional derivative is $$u^\,_{\{1,y\}}= \text{ grad }u\cdot \{1,y\}= \{u_x,u_y\}\cdot \{1,y\} = u_x+yu_y;$$

  • Since characteristic lines are defined by the equation $u(x,y) = C,$ then the full differential $$\text{ d}u = u_x \text{ d}x + u_y \text{ d}y = 0$$ on this lines is zero, and $$\frac{\text{ d}y}{\text{ d}x} = -\dfrac{u_x}{u_y}.$$ Taking in account PDE $(OP.4),$ this leads to the ODE $(OP.5).$

  • To get the solution from $(OP.7),$ Maclaurin series can be used. However, the kernel of the representation $(5)$ continuously depends of $\lambda$ and looks more suitable if $\lambda$ is not integer.