How to intuitively arrive at the total derivative limit and the jacobian matrix? I'm following this PDF and I need to understand how to arrive at the definition of total derivative geometrically.
For now, what I understand is that, from the original definition of derivative:
$$\lim_{h\to 0} \frac{f(p+h)-f(x)}{h} = m \implies \lim_{h\to 0}\frac{f(x+p)-f(x)-mh}{h} = 0$$
Which is the same as
$$\lim_{|h|\to 0}\frac{|f(x+p)-f(x)-mh|}{|h|} = 0$$
and the advantage of this is that we can know define the existence of the derivative, but now for a multivariable function $f$ from $\mathbb R^m \to \mathbb R^n$. Besides not giving any geometrial interpretation, I quite accept this result.
But then, the PDF starts just defining what the total derivative would be. I need an intuitive way to arrive at this, and then, a way to start from the limit definition of the total derivative and arrive at the jacobian matrix.

Also, how it was historically derivated for the first time?

 A: For functions $f:\>{\mathbb R}\to{\mathbb R}$, ${\mathbb R}\to{\mathbb R}^n$ or ${\mathbb C}\to{\mathbb C}$ we can define the derivative in a straightforward way as a limit
$$f'(x):=\lim_{h\to 0}{f(x+h)-f(x)\over h}\ ,\tag{1}$$
because the increment variable $h$ takes values in a field (${\mathbb R}$ or ${\mathbb C}$) and can be put into the denominator of the expression on the right hand side of $(1)$.
In a multivariate situation this is no longer possible: We cannot divide by vectors. Therefore we have to explain the derivative in a "denominator-free" way. The original setup $(1)$ can be reformulated as follows: Given $f$ and a point $x$ in the domain of $f$ there is a constant $A$ such that
$$f(x+h)-f(x)=A\>h+ o\bigl(|h|\bigr)\qquad(h\to 0)\ .\tag{2}$$
This number $A$ is called the derivative of $f$ at $x$, and is denoted by $f'(x)$. The formula $(2)$ can be interpreted intuitively as follows: For "small" increments $h$ "attached" at $x$ the resulting increment $\Delta f=f(x+h)-f(x)$ of the function value is in first approximation a constant multiple of $h$. This principle is capable of being transferred to a multivariate situation: For "small" increment vectors $h$ attached at $x$  the resulting increment vector $\Delta f=f(x+h)-f(x)$ is in first approximation a linear function of the increment vector $h$. In vector notation this looks as follows:
$${\bf f}({\bf x}+{\bf h})-{\bf f}({\bf x})=A{\bf h}+o\bigl(|{\bf h}|\bigr)\qquad({\bf h}\to{\bf0})\ .\tag{3}$$
The linear map $A:{\mathbb R}^n\to{\mathbb R}^m$ appearing here is called the derivative (not the "total derivative") of ${\bf f}$ at ${\bf x}$, and is denoted by $d{\bf f}({\bf x})$, or similar. The equation $(3)$ can then be written in the form
$${\bf f}({\bf x}+{\bf h})-{\bf f}({\bf x})=d{\bf f}({\bf x}).{\bf h}+o\bigl(|{\bf h}|\bigr)\qquad({\bf h}\to{\bf0})\ .\tag{4}$$
When $f$ happens to be a real-valued function the formula $(4)$ assumes the form
$$f({\bf x}+{\bf h})-f({\bf x})=\nabla f({\bf x})\cdot{\bf h}+o\bigl(|{\bf h}|\bigr)\qquad({\bf h}\to{\bf0})\ .$$
A: I would start with the gradient.
Consider $\phi:\mathbb{R}^{N}\rightarrow\mathbb{R}$
$$\delta \phi=\phi(x_{1}+\delta x_{1},\cdots,x_{N}+\delta x_{N})-\phi(x_{1},\cdots,x_{N})$$
Taking a first order Taylor expansion
$$\delta \phi = \frac{\partial \phi}{\partial x_{1}}\delta x_{1}+\cdots+\frac{\partial \phi}{\partial x_{N}}\delta x_{N}+\mathcal{O}(\delta x^{2})$$
Recognizing the above as an inner product  
$$\delta \phi  = \nabla\phi\cdot\langle \delta x_{1},\cdots \delta x_{N} \rangle $$
Now I divide by the length between the two points in the domain
$$\frac{\delta \phi}{\delta x}=\nabla\phi\cdot\langle \delta x_{1},\cdots \delta x_{N} \rangle\frac{1}{\delta x} $$
Where $\delta x = \sqrt{\sum_{j=1}^{N}\delta x_{j}^{2}}$
We then define
$$\hat{n}=\langle \delta x_{1},\cdots \delta x_{N} \rangle\frac{1}{\delta x}$$
$$\frac{\delta \phi}{\delta x}=\nabla\phi \cdot \hat{n}$$
Where you would now like to think of $\nabla\phi$ as the derivative, it carries all of the information of all directional derivatives and you simply need to take an inner product with any normal vector to get a directional derivative.
You can run this same argument for a function that has an $N$ dimensional domain and an $M$ dimensional codomain. There is a similar interpretation for the Jacobian. 
A: I guess it will be hard to answer this without knowing what you mean by "geometrically".  How geometric is a linear map between arbitrary vector spaces to you?
The intuition for the derivative is that it is the best linear approximation to the function at a point.
So, for me,  $f(\vec{x}+\vec{h}) \approx f(\vec{x})+ Df\big|_\vec{x}(\vec{h})$,  where $Df\big|_\vec{x}$ is a linear function, is enough intuition.
To make this intuition rigorous, we need to know in what sense the approximation holds.  We know it should get better and better as $\vec{h}$ gets smaller.  The error between the approximation and the actual function value should be much smaller than the difference in the inputs.  One way to formulate such a condition is to say
$$f(\vec{x}+\vec{h}) = f(\vec{x})+Df\big|_\vec{x}(\vec{h})+\textrm{Error}(\vec{x},\vec{h})$$
where $Df\big|_\vec{x}(\vec{h})$ is linear and $\displaystyle \lim_{\vec{h} \to 0} \frac{\left| \textrm{Error}(\vec{x},\vec{h}) \right| }{\left| \vec{h} \right|} = 0$.  This is equivalent to your definition.
Now lets see if we can see the connection to the jacobian matrix.  The columns of the matrix of $Df\big|_\vec{x}$ should be what I get from plugging the standard basis vectors into the approximation above.  Problem is, those standard basis vectors have big norms.  So just plug in scalar multiples of them instead, and see what happens.  In other words, I use a displacement of $h\vec{e_j}$ where $h$ is a scalar, and $e_j$ is a basis vector.
$\begin{align*}
f(x+he_j) &= f(x)+Df\big|_x(he_j) + \textrm{Error}(x,he_j)\\
f(x+he_j) &= f(x)+hDf\big|_x(e_j) + \textrm{Error}(x,he_j)\\
\frac{f(x+he_j)-f(x)}{h} &= Df\big|_x(e_j) + \frac{\textrm{Error}(x,he_j)}{h}\\
\lim_{h \to 0 }\frac{f(x+he_j)-f(x)}{h} &= Df\big|_x(e_j)
\end{align*}$
Which gives you that the entries in the matrix representing $Df\big|_x$ are the appropriate partial derivatives.
