Recall the theorems from calculus 3:

Theorem 2: If $g(t)=(x(t),y(t))$ is a smooth path, $z=f(x,y)$ is differentiable with nonzero gradient and the velocity $g'(t)$ is orthogonal with the gradient $\nabla f $ at $g(t)$, then $f(g(t))$ is constant, that is, the range of $g$ lies in a level curve of $f$.

We can use these theorems to build approximations to the level and gradient curves of $f$ through the point $(x_0,y_0)$. It's important that the gradient of $f$ be nonzero in a neighborhood of $(x_0,y_0)$. For then we can build the ascent (descent) curve of $f$ starting from $(x_0,y_0)$, by choosing a small positive step size $dt=1/N$, and proceed by stepping $(x_{n+1},y_{n+1})=(x_n,y_n)+(-)\dfrac{\nabla f(x_n,y_n)}{|\nabla f(x_n,y_n)|}\,dt$ for $n=0,1,2, \cdots, N$. Similarly, for the level curve (just rotating the gradient vector 90 degrees and then stepping). This method (which is just the basic Euler method for numerical solution of a differential equation) works pretty well to obtain a local image of the level and gradient curves of $f$. By choosing a list of starting points $(x_0,y_0)$, we can piece together these patches to get a bigger picture. Try the Sagecell interact below to experiment with this. If you choose the 3d option t3d, give the image time to form before manipulating it.