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}}. $