Update
Thanks to the comment below. I fixed a silly error I had. But now after fixing the error, I found I am not able to solve this PDE.
So the question really becomes, can this pde be solved analytically? This explains why Maple and Mathematica could not solve it.
Here is the corrected steps I did, and showing where I stuck.
Solve
$$
xu_{x}+yu_{y}=-4xyu
$$
with Initial conditions
$$
u\left( x,0\right) =e^{-x^{2}}
$$
Using the method of characteristics (with parameterization), the systems of
characteristic lines are (from the PDE itself)
\begin{align}
\frac{dx}{ds} & =x\tag{1}\\
\frac{dy}{ds} & =y\tag{2}\\
\frac{du}{ds} & =-4xyu\tag{3}
\end{align}
With initial conditions
\begin{align*}
x\left( 0\right) & =\xi\\
y\left( 0\right) & =0\\
u\left( 0\right) & =e^{-\xi^{2}}
\end{align*}
Solving (1) gives
$$
x=C_{1}e^{s}
$$
At $s=0,x=\xi$. Hence $C_{1}=\xi$ and
$$
x=\xi e^{s}
$$
Solving (2) gives
$$
y=C_{1}e^{s}
$$
At $s=0,y=0$. Hence $C_{1}=0$ and
$$
y=0
$$
From the above we see that
$$
\xi=xe^{-s}
$$
But no second equation exist to solve for $s$, this is because $y=0$.
Equation (3) becomes
$$
\frac{du}{ds}=-4xyu
$$
But $y=0$, hence the above becomes
\begin{align*}
\frac{du}{ds} & =0\\
u & =u\left( 0\right) \\
& =e^{-\xi^{2}}
\end{align*}
But can not convert back to $x,y$ since $\xi=xe^{-s}$ and do not have an
equation to solve for $s$. Hence no solution.
Original question (has error) kept for reference
I solved this pde by hand
$$
x u_x + y u_y = -4 x y u
$$
With initial conditions $u(x,0)=e^{-x^2}$. The strange thing is that the solution I get verifies the PDE OK. But when $y=0$ the solution does not give $e^{-x^2}$ which is the initial conditions.
I am not able to see where my mistake is. Went over it few times.
I'll show below the steps I used to solve the above. Could someone show where the error is? I am assuming the solution is wrong (even though it satisfies the PDE) because it should also satisfy the initial conditions.
The computer can not solve this PDE for some reason. (Tried Maple and Mathematica).
Using the method of characteristics (with parameterization), the systems of
characteristic lines are (from the PDE itself)
\begin{align}
\frac{dx}{ds} & =x\tag{1}\\
\frac{dy}{ds} & =y\tag{2}\\
\frac{du}{ds} & =-4xyu\tag{3}
\end{align}
With initial conditions (at $s=0$)
\begin{align*}
x\left( 0\right) & =\xi\\
y\left( 0\right) & =0\\
u\left( 0\right) & =e^{-\xi^{2}}
\end{align*}
Solving (1) gives
\begin{align*}
\ln x & =s+x\left( 0\right) \\
& =s+\xi\\
x & =\xi e^{s}
\end{align*}
Solving (2) gives
\begin{align*}
\ln y & =s+y\left( 0\right) \\
& =s\\
y & =e^{s}
\end{align*}
From the above we see that
\begin{align}
\xi & =\frac{x}{y}\tag{4}\\
s & =\ln y\nonumber
\end{align}
Equation (3) becomes
$$
\frac{du}{u}=-4xyds
$$
But from earlier we found that $x=\xi e^{s}$ and $y=e^{s}$, hence the above becomes
\begin{align*}
\frac{du}{u} & =-4\left( \xi e^{s}\right) \left( e^{s}\right) ds\\
& =-4\xi e^{2s}ds
\end{align*}
Integrating gives
$$
\ln u=-2\xi e^{2s}+u\left( 0\right)
$$
Therefore
\begin{align*}
u & =u\left( 0\right) e^{-2\xi e^{2s}}\\
& =e^{-\xi^{2}}e^{-2\xi e^{2s}}
\end{align*}
Changing back to $x,y$ using (4) gives
\begin{align*}
u\left( x,y\right) & =e^{-\left( \frac{x}{y}\right) ^{2}}e^{-2\frac
{x}{y}e^{2\ln y}}\\
& =e^{-\left( \frac{x}{y}\right) ^{2}}e^{-2xy}
\end{align*}
Verified the solution OK, but there is something wrong. At $y=0$ I should get
$u\left( x,0\right) =e^{-x^{2}}$. But this above gives $u\left( x,0\right)
=0$.
Best Answer
You get, first disregarding the initial conditions, in a non-parametric calculation, $y=c_1x$ and $$ \frac{dx}{x}=\frac{dy}{y}=\frac{d(xy)}{2xy}=\frac{du}{-4xyu} $$ so that $ue^{2xy}=c_2=\phi(c_1)=\phi(\frac yx)$.
From the first identity you see that the line $y=0$ is already a characteristic curve, the second identity then demands that $u$ be constant to $\phi(0)$ on that line. So you can not give a non-constant initial condition on that line.
You made no error, this task has no solution.
If the initial condition were, as an example, modified to $u(x,1)=e^{-x^2}$, one would get $\phi(\frac1x)=e^{-x^2+2x}$ so that then $$ u(x,y)=\exp\left(-\frac{x^2}{y^2}+\frac{2x}y-2xy\right). $$