This is a community-wiki answer trying to remove this question from the unanswered queue.
Find $\nabla f$ at the point $P_0 = (0, 1)$.
For the implicitly defined $z = f(x,y)$ by
$$
x^3 +y^3 +z^3 + xyz =9,\tag{1}
$$
the general procedure from a standard multivariable calculus course is taking partial derivatives w.r.t. $x$ and $y$ in (1) respectively, treating $z$ as a function of $x$ and $y$ (chain rule kicks in). Therefore we have:
$$
3x^2 + 3z^2 \frac{\partial z}{\partial x} + yz + xy\frac{\partial z}{\partial x} = 0,
\\
3y^2 + 3z^2 \frac{\partial z}{\partial y} + xz + xy\frac{\partial z}{\partial y} = 0.
$$
Solving for $\partial z/\partial x$ and $\partial z/\partial y$ we have:
$$
\frac{\partial z}{\partial x} = -\frac{3x^2+yz}{3z^2+xy},\quad\text{and}\quad
\frac{\partial z}{\partial y} = -\frac{3y^2+xz}{3z^2+xy}.
$$
Plugging $x= 0,y=1$ and $z=2$:
$$
\nabla f = \left(\frac{\partial z}{\partial x} , \frac{\partial z}{\partial y} \right)
= \left(-\frac{1}{6} , -\frac{1}{4} \right).
$$
Find the rate and direction of the steepest increase of $f$ at $P_0$.
At a point of a surface $z = f(x,y)$, the steepest increase direction coincides with the direction of the gradient. Rate is the magnitude of the gradient: $ |\nabla f |$.
Find $(D_u f)(P_0)$, where $u = (\mathbf{i}+\mathbf{j})/\sqrt{2}$.
Using the results above:
$$
(D_u f)(P_0) = \nabla f\cdot u\big|_{P_0} = \left(-\frac{1}{6} , -\frac{1}{4} \right) \cdot \left( \frac{1}{\sqrt{2}} , \frac{1}{\sqrt{2}} \right) = -\frac{5}{12\sqrt{2}}.
$$