In general, if a function from [itex]R^m[/itex] to [itex]R^n[/itex] is "differentiable at [itex]\vec{v_0}[/itex]" if and only if there exist a linear function, L, from [itex]R^m[/itex] to [itex]R^n[/itex] and a function [itex]\epsilon[/itex] from [itex]R^m[/itex] to [itex]R^n[/itex] such that [itex]f(\vec{v})= f(\vec{v_0})+ L(\vec{v}- \vec{v_0})+ \epsilon(\vec{v})[/itex] and [itex]\lim_{\vec{v}\to \vec{0}} \epsilon(\vec{v}/|\vec{v}|[/itex][itex]= 0[/itex]. In this case, we say that L is the "derivative of f with respect to [itex]\vec{v}[/itex] at [itex]vec{v_0}[/itex]".
Notice that L is a linear operator, not a number or even a vector. In the case f from [itex]R^1[/itex] to [itex]R^1[/itex], a real valued function of a single real variable as in Calculus I, a linear function from [itex]R^m[/itex] to [itex]R^n[/itex] is a multiple: [itex]mx[/itex] for some number m. In calculus I, we think of that number as being the derivative.
In the case of f from [itex]R^m[/itex] to [itex]R^1[/itex], a real valued function of several real variables, we can represent a linear function from [itex]R^m[/itex] to [itex]R^1[/itex] as a "dot product"- any linear function L(<x, y, z>) can be written as ax+ by+ cz= <a, b, c>[itex]\cdot[/itex]<x, y, z>. In that case, we can identify the derivative with that vector- which is, in fact, the gradient.
That is the case in your example. f is a real valued function of a variable in [itex]R^3[/itex]. Strictly speaking, its derivative is the linear transformation corresponding to taking the dot product of the position vector [itex]x\vec{i}+ y\vec{j}+ z\vec{k}[/itex] with the vector
[tex]\frac{\partial f}{\partial x}\vec{i}+ \frac{\partial f}{\partial y}\vec{j}+ \frac{\partial f}{\partial z}\vec{k}[/tex]
which is, as you say, the gradient of f.
(Notice that we need "vector space" structure in the range space of f because we need to be able to subtract f(x+h)- f(x). We don't need to do "arithmetic" in the domain space because all we do there is take the limit as |h| goes to 0. That's why we tend to talk about "functions of several variables" rather than "functions of a vector variable".)
In the general case of f from [itex]R^m[/itex] to [itex]R^n[/itex] with neither m nor n equal to 1, a linear function from [itex]R^m[/itex] to [itex]R^n[/itex] can be written as a matrix multiplication and we can identify the derivative with that matrix. In the case of [itex]\vec{f}(\vec{v})= f_x(x, y, z)\vec{i}+ f_y(x,y,z)\vec{j}+ f_z(x,y,z)\vec{k}[/itex] (where \vec{v}= x\vec{i}+ y\vec{j}+ z\vec{k}[/itex]) that matrix is
[tex]\begin{bmatrix}\frac{\partial^2 f_x}{\partial x^2} & \frac{\partial^2 f_y}{\partial x\partial y} & \frac{\partial^2 f_z}{\partial x\partial z} \\ \frac{\partial^2 f_x}{\partial y\partial x} & \frac{\partial^2 f_y}{\partial y^2} & \frac{\partial^2 f_z}{\partial y\partial z}\\ \frac{\partial^2 f_z}{\partial z\partial x} & \frac{\partial^2 f_z}{\partial z\partial y} & \frac{\partial^2 f_z}{\partial z^2}\end{bmatrix}[/tex]