I also do not know how to tell Mathematica to give step by step. May be the Wolfram Alpha PRO version would show the steps? But this is how I would solve this by hand.
Solve, where $u\equiv u\left( x,y\right) $
$$
\frac{3x^{2}}{\frac{\partial u}{\partial x}}+\frac{3y^{2}}{\frac{\partial
u}{\partial y}}=-1
$$
We assume that
\begin{equation}
u=F_{1}\left( x\right) +F_{2}\left( y\right) \tag{1}
\end{equation}
Plugging the above into the PDE gives
\begin{align*}
\frac{3x^{2}}{F_{1}^{\prime}}+\frac{3y^{2}}{F_{2}^{\prime}} & =-1\\
\frac{3x^{2}}{F_{1}^{\prime}} & =-1-\frac{3y^{2}}{F_{2}^{\prime}}
\end{align*}
Where $F_{1}^{\prime}\equiv\frac{dF_{1}}{dx},F_{2}^{\prime}=\frac{dF_{2}}{dy}
$. Since each side depends on different variables and both are equal, then
they must be both constant. Say $c$. We obtain two differential equations
\begin{align*}
\frac{3x^{2}}{F_{1}^{\prime}} & =c\\
-1-\frac{3y^{2}}{F_{2}^{\prime}} & =c
\end{align*}
Hence
\begin{align*}
F_{1}^{\prime} & =\frac{3x^{2}}{c}\\
F_{2}^{\prime} & =\frac{-3y^{2}}{1+c}%
\end{align*}
The solution to the first ODE is
$$
F_{1}\left( x\right) =\frac{x^{3}}{c}+c_{1}
$$
The solution to the second ODE is
$$
F_{2}\left( y\right) =\frac{-y^{3}}{1+c}+c_{2}
$$
Hence from (1)
\begin{align*}
u & =F_{1}\left( x\right) +F_{2}\left( y\right) \\
& =\frac{x^{3}}{c}+c_{1}-\frac{y^{3}}{1+c}-c_{2}
\end{align*}
Let $c_{1}-c_{2}=\tilde{c}_{2}$, therefore the final solution is
$$
u(x,y) =\frac{x^{3}}{c}-\frac{y^{3}}{1+c}+\tilde{c}_{2}
$$
To verify, plugin the above solution into the PDE
$$
\frac{3x^{2}}{\frac{\partial u}{\partial x}}+\frac{3y^{2}}{\frac{\partial
u}{\partial y}}=\frac{3x^{2}}{\frac{3x^{2}}{c}}+\frac{3y^{2}}{\frac{-3y^{2}
}{1+c}}=c-\left( 1+c\right) =-1
$$
QED
update
To answer comment asking for Mathematica code. Here is a possible way. But I am sure this can be done better. Write down the pde and do the replacement
ClearAll[u, x, y, f1, f2, c]
eq = 3 x^2/D[u[x, y], x] + 3 y^2/D[u[x, y], y] == -1;
eq2 = eq /. u -> Function[{x, y}, f1[x] + f2[y]]
Now do the separation of variables. Had to do this by looking at the answer above. Automating this step is hard for me, in other words, it is left as an exercise.
sol1 = f1[x] /. First@DSolve[eq2[[1, 1]] == c, f1[x], x];
sol2 = f2[y] /. First@DSolve[-eq2[[1, 2]] - 1 == c, f2[y], y];
mySolution = sol1 + sol2
Now just plug the above solution back into the PDE to verify
eq /. u -> Function[{x, y}, Evaluate@mySolution]
QED
