This is in the context of explaining tangent lines as linear approximations, where the function of the tangent line is given by $T(x) = f(x_0) + f'(x_0)(x-x_0)$. The author mentions that $\left\lvert \frac{f(x) - f(x_0)}{x - x_0} - f'(x_0) \right\lvert \to 0$ as $x \to x_0$ implies $\left\lvert \frac{f(x) - T(x)}{x -x_0} \right\lvert \to 0$ as $x \to x_0$, which I have no problem with: that's simply a linear approximation which minimizes the distance between the approximation $T(x)$ and the function $f(x)$.
What I can't figure out is why the author says it implies$f(x_0 + h) - f(x_0) = f'(x_0)h + \epsilon h$, where $h = x - x_0$ and $\epsilon$ goes to $0$ as $h$ goes to $0$, which is meant to be the beginning of a generalization for multivariable cases.
I would actually have expected something like $f(x_0 + h) - f(x_0) = f'(x_0)h + \epsilon(h)$, where $\epsilon(h)$ is $|f(x_0 +h) - T(x_0+h)|$, i.e., the error function describing the distance between the approximation and the function approximated. But instead of $\epsilon(h)$, the author uses $\epsilon h$. A possibly important fact here is that $\lim_{h \to 0} f(x) = L$ implies that $f(x) = L - g(x)$ where $g(x)$ is some function s.t. $\lim \limits_{h \to 0} g(x) = 0$, and this probably has something to do for why there is an $\epsilon h$ there, but I can't understand why exactly the author chose $\epsilon h$ (it seems to me there are many possible choices for $g(x)$) or how exactly he gets it in the equation: the author explicitly says it is nothing more than a manipulation of $\left\lvert \frac{f(x) - T(x)}{x -x_0} \right\lvert \to 0$ as $x \to x_0$.
Any clues?
The problem is from page 709 of Elementary Real Analysis by Thomson, Bruckner and Bruckner, and they use $$f(x_0 + h, y_0 + k) - f(x_0, y_0) = f_1(x_0, y_0)h + f_2(x_0, y_0)k + \epsilon(h,k)\sqrt{h^2 + k^2}$$ in the case of a two variable function.
I think writing $f(x_0 + h) - f(x_0) = f'(x_0)h + \epsilon h$ is misleading because it suggests that $\epsilon$ is some constant. But in fact you have $$f(x_0 + h) - f(x_0) = f'(x_0)h + \epsilon(h) h \tag{1}$$ with $$\epsilon(h) = \frac{f(x_0 +h) - f(x_0)}{h} - f'(x_0) , \tag{2}$$ i.e. $\epsilon$ is a function defined on $\{ h \in \mathbb R \mid x_0 + h \text{ lies in the domain of } f \}$. Formula $(1)$ is absolutely trivial - so what is its benefit? The point is that $$T(x) = f(x_0) + f'(x_0)(x-x_0) = f'(x_0)x + (f(x_0) - f'(x_0)x_0)$$ is a linear function which is the best linear approximation of $f$ in the sense that the relative error $$\epsilon(h) = \frac{f(x_0 +h) - f(x_0)}{h} - f'(x_0) = \frac{f(x_0 +h) - T(x_0+h)}{h}$$ goes to $0$ as $h$ goes to $0$. No other linear $\tilde T(x) = ax + b$ has this property.