How does Mathematica solve a certain differential equation?

@Nasser's method seems to be what Mathematica does internally. The following gives hints of the steps taken, and one can see the elements of Nasser's solution. Some is left to the user to guess. The first code modifies Integrate to print itself out; the second uses Trace to see the calls to Integrate.

Block[{DSolve`print = Print},                                (* internal hook *)
 Internal`InheritedBlock[{Integrate}, Unprotect[Integrate];  (* alter Integrate to print itself *)
  i : Integrate[___] /; ! TrueQ[$in] := 
   Block[{$in = True}, Print["***Integrate= ", HoldForm[i]]; i];
  Protect[Integrate];
  DSolve[3*x^2/D[u[x, y], x] + 3*y^2/D[u[x, y], y] == -1, u, {x, y}]
  ]]

This does something similar, but the Integrate commands are output at the end, not at the time they are called:

Block[{DSolve`print = Print},
 Trace[
  DSolve[3*x^2/D[u[x, y], x] + 3*y^2/D[u[x, y], y] == -1, u, {x, y}],
  _Integrate,
  TraceInternal -> True
  ]]

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