$\textbf{Comment}$: For all that is about to follow, I have better notes on here derivative as a linear map (pg 4).
I can start you off with $f: \mathbb{R}^n \to \mathbb{R}^m$ where $m=1$ which is the directional derivative. You can generalize the argument from this special case to the standard derivative.
$$f \in C^1(U \subset \mathbb{R}^n,\mathbb{R}) \Rightarrow Df(p) \in L( \mathbb{R}^n, \mathbb{R})$$
i.e $Df(p)$ is a linear map from $\mathbb{R}^n$ to $\mathbb{R}$. Given $(x_1,...,x_n)$ are the coordinate functions on $U$ and $\alpha:(- \epsilon, \epsilon) \to U$ is a differentiable curve (with $\alpha(0) = p$) then it is natural to have,
\begin{align} Df(p) \cdot \alpha'(0)&= \frac{d}{dt} \Bigr|_{t = 0} (f \circ \alpha)(t) \\ \\ & = \nabla f(p) \cdot \alpha'(0) \\ \\ & \Rightarrow Df(p) = \nabla f(p) \end{align}
The reason why this is natural to want is that $f$ is differentiable at $p$ i.e by the inverse function theorem, it is a local diffeomorphism i.e it maps curves to curves.
Hence, the curve $\alpha$ gets mapped to the curve $f \circ \alpha$ and so you would want the range vector to $\alpha$ at $0$ to be the tangent vector to $f \circ \alpha$ at $0$. We now recall the definition of differentiability. We say $g: U \subset \mathbb{R}^n \to V \subset \mathbb{R}^m$ is differentiable at $p \in U$ if there exists a linear map $\lambda: \mathbb{R}^n \to \mathbb{R}^m$ such that the following holds,
$$ \lim_{x \to p} \frac {\|f(x) - f(p) - \lambda(x-p)\|_m}{\|x-p\|_n} = 0$$
You can check that such a $\lambda$ must be unique. Here $x-p$ is just a tangent vector to $p$. Hence, you can just check that,
$$ \nabla f(p) \cdot \alpha'(0) = \sum_{j=1}^{n} \frac{\partial f}{\partial x_j}(p) \cdot \alpha_j'(0), \textbf{where} \ \ \alpha(t) = \begin{pmatrix} \alpha_1(t) \\ \vdots \\ \alpha_n(t) \end{pmatrix}$$
is $\lambda$. This is shown in any multivariable calculus book as the linearization of a function. They will typically write,
\begin{align} z(x,y)& \approx f_x(a,b)(x-a) + f_y(a,b)(y-b) + f(a,b) \\ \\ & = \nabla f(a,b) \cdot \begin{pmatrix} x-a \\ y-b \end{pmatrix} + f(a,b)\end{align}
This argument is done for $f: \mathbb{R}^2 \to \mathbb{R}$. Here $z(x,y)$ is an approximation of $f(x,y)$ for $(x,y)$ near $(a,b)$. Therefore we have the following,
$$f(x,y) - f(a,b) - \nabla f(a,b) \cdot\begin{pmatrix} x-a \\ y-b \end{pmatrix} = \epsilon(x,y)$$
If you accept this special case to be true, you see the extension almost immediately. This time you get a matrix for $Df(p)$. Well you had one last time, it was just $1 \times n$ or $n \times 1$. So similarly you guess that $Df(p) = \nabla f(p)$ as above, and now when checking, you'll see it amounts to this special case in each entry. Applying the triangle inequality (inductively) you get the result.