Step 1: The given PDE \(u_x + y u_y = 0\) is a linear first order PDE of the form \(Pu_x + Qu_y = R\) with \(P=1\), \(Q=y\), \(R=0\). By Lagrange's method, the auxiliary equations are
\[\frac{dx}{1}=\frac{dy}{y}=\frac{du}{0}\]Step 2: From \(du = 0\), we get \(u = c_1\), a constant along each characteristic curve.
Step 3: From \(\frac{dx}{1}=\frac{dy}{y}\), separate variables: \(\frac{dy}{y} = dx\). Integrating gives \(\ln y = x + k\), so
\[y e^{-x} = c_2\]Step 4: The general solution of the PDE is \(u = F(c_2) = F(y e^{-x})\) for an arbitrary function \(F\).
Step 5: Apply the initial condition \(u(0,y) = y^3\). At \(x=0\), \(c_2 = y\), so \(F(y) = y^3\), which means \(F(s) = s^3\) for any argument \(s\).
Step 6: Substitute \(s = y e^{-x}\) back into \(F\):
\[u(x,y) = (y e^{-x})^3 = y^3 e^{-3x}\] \[\boxed{u(x,y) = y^3 e^{-3x}}\]