I am trying to solve the PDE
\begin{align} xu_x + yu_y = 2xy \end{align}
using the method of characteristics. So the characteristic equations are
\begin{align} \frac{dx}{x} = \frac{dy}{y} = \frac{du}{2xy} \end{align}
Considering the first and second equation we get
\begin{align} x = c_1y,\quad c_1\,\text{constant} \end{align}
Considering the second and third equation we get
\begin{align} 2xy = u + c_2,\quad c_2\text{ constant} \end{align}
therefore the general solution is given by \begin{align} u(x, y) = 2xy - f\left(\frac{x}{y}\right), \quad f\text{ arbitrary function} \end{align}
However if I substituite this solution in the PDE I get \begin{align} xu_x + yu_y = 4xy\,(?) \end{align}
By the other hand, if I write down the characteristic equations parametrically, i.e
\begin{align} \frac{dx}{ds} &= x \\ \frac{dy}{ds} &= y \\ \frac{du}{ds} &= 2xy \end{align}
from the first two I get \begin{align} x = Ae^s \quad\text{ and }\quad y = Be^s,\quad A,\,B \text{ constants} \end{align} Eliminating the parameter $s$ we get \begin{align} (*)\qquad \frac{x}{y} = \frac{A}{B} = C,\quad C\text{ constant} \end{align} Then to solve for $u$ \begin{align} \frac{du}{ds} &= 2xy = 2ABe^{2s}\quad \text{ since } x = Ae^s,\quad y = Be^s \end{align} Integrating \begin{align} \quad u &= e^{2s}AB + D,\quad D \text{ constant} \\ \end{align} Now since $xy = ABe^{2s}$, and using (*) we put $D = f(x/y)$, it follows \begin{align} u = xy + f\left(\frac{x}{y}\right) \end{align} and this is the solution since subtituing in the PDE I get an identity.
Why am I getting the incorrect answer with the first method?
$$xu_x + yu_y = 2xy$$ You wrote :
$\frac{dx}{x} = \frac{dy}{y} = \frac{du}{2xy}\quad $is OK.
$x = c_1y\quad $ is OK.
$2xy = u + c_2\quad$ is false.
From $\frac{dy}{y} = \frac{du}{2xy}\quad\implies\quad 2x\,dy=du$.
Then you try to integrate. But $\quad \int 2x\,dy\neq 2xy+c\quad$ because $x$ is not constant. Without writting explicitly the function $x(y)$ you cannot integrate. ,
HINT : $\quad 2c_1y\,dy=du$