Directional Derivative in $\mathbb{R}^2$ 
Definition Let $f: \mathbb{R}^n \to \mathbb{R}$ and $u \in \mathbb{R}^n$ be a unit vector. The directional derivative of $f$ in the direction of $u$ is

$$D_uf(j) = \displaystyle\lim_{t \to 0} \frac{f(j+tu) - f(j)}{t}$$
provided that this limit exists.

I am preparing for my summer exit exams in my grad program. I saw the following problem in my notes, but with no solution. So I attempted it myself. I am looking for solution verification. Please correct me if I have made any mistakes along the way.

Example Suppose $f(x,y) = x^2+3xy+4y^2$ and $j =(2,1)$ and $u = \langle \frac{3}{5}, - \frac{4}{5} \rangle$. To find $D_uf(j)$, we find the limit.
$$D_uf(j) = \displaystyle\lim_{t \to 0} \frac{f(j+tu) - f(t)}{t} = \displaystyle\lim_{t \to 0} \frac{f((2,1)+ t(\frac{3}{5}, \frac{-4}{5}))-f(2,1)}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{f(2+\frac{3}{5}t, 1 - \frac{4}{5}t)-f(2,1)}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{(2+\frac{3}{5}t)^2+3(2+\frac{3}{5}t)(1-\frac{4}{5}t) +4(1-\frac{4}{5}t)^2-[2^2-3(2)(1)+4(1)^2]}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{( \frac{1337}{100}t^2+\frac{23}{5}t+14) - 2}{t} = \displaystyle\lim_{t \to 0} \frac{\frac{1337}{100}t^2+\frac{23}{5}t+12}{t}$$
$$= \displaystyle\lim_{t \to 0} 13.37t + \displaystyle\lim_{t \to 0} \frac{23}{5} + \displaystyle\lim_{t \to 0} \frac{12}{t} = \frac{23}{5} + \infty$$
How do I get rid of this pesky $\frac{12}{t}$? Perhaps I have made an error somewhere.

EDIT / UPDATE:
This only works if we have $f(x,y) = x^2-3xy+4y^2$, where the second term is negative instead of positive. Doing so yields
$$D_uf(j) = \displaystyle\lim_{t \to 0} \frac{f(j+tu) - f(t)}{t} = \displaystyle\lim_{t \to 0} \frac{f((2,1)+ t(\frac{3}{5},\frac{-4}{5}))-f(2,1)}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{f(2+\frac{3}{5}t, 1 - \frac{4}{5}t)-f(2,1)}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{(2+\frac{3}{5}t)^2-3(2+\frac{3}{5}t)(1-\frac{4}{5}t) +4(1-\frac{4}{5}t)^2-[2^2-3(2)(1)+4(1)^2]}{t}$$
$$= \displaystyle\lim_{t \to 0} \frac{\frac{109}{25}t^2 -t +2 - [2]}{t} = \displaystyle\lim_{t \to 0} = \displaystyle\lim_{t \to 0} \Big( \frac{109}{25}t - 1 \Big) = -1$$
 A: Maybe you can recognize that
$Df(\mathbf{x}_0)[\mathbf{u}] 
=g'(0)$ with the scalar-valued function
$g(t) = f(\mathbf{x}_0 + t\mathbf{u} )$.
Denote $\mathbf{x}= \mathbf{x}_0 + t\mathbf{u}$.
It is simple to show using chain rule that
$g'(t) = \nabla_\mathbf{x}f(\mathbf{x}_0 + t\mathbf{u}):\mathbf{u}$
from which you can deduce
$$
Df(\mathbf{x}_0)[\mathbf{u}]= g'(0)
= \nabla_\mathbf{x}f(\mathbf{x}_0):\mathbf{u}
$$
In your application
$$
\nabla_\mathbf{x}f(\mathbf{x}_0)
=
\begin{pmatrix}
2x_0+3y_0 \\
8y_0+3x_0
\end{pmatrix}
=
\begin{pmatrix}
7 \\
14
\end{pmatrix}
$$
The directional derivative is
$\frac15 (7\cdot 3-14\cdot 4)=-7$
Here the dot colon indicates the inner product between vectors.
A: Observe that the proposed function $f:\mathbb{R}^{2}\to\mathbb{R}$ is a polynomial, hence $C^{\infty}$.
In particular, it is differentiable.
We can then state the following auxiliary result:
Proposition
If the function $f:\mathbb{R}^{2}\to\mathbb{R}$ is differentiable, then its derivative along the direction $v\in\mathbb{R}^{2}$ at the point $p\in\mathbb{R}^{2}$ is given by the corresponding formula:
\begin{align*}
D_{v}f(p) = f'(p)v
\end{align*}
Proof
Due to the fact that $f$ is differentiable, the next limit exists:
\begin{align*}
\lim_{h\to 0}\frac{|f(p + h) - f(p) - f'(p)h|}{\|h\|} = 0
\end{align*}
Hence, if we make the change of variable $h = tv$, it results that
\begin{align*}
\lim_{t\to 0}\frac{1}{\|v\|}\frac{|f(p + tv) - f(p) - f'(p)tv|}{|t|} = \lim_{t\to 0}\frac{1}{\|v\|}\left|\frac{f(p + tv) - f(p)}{t} - f'(p)v\right| = 0
\end{align*}
whence we conclude that $D_{v}f(p) = f'(p)v$. In particular, since $f'(p)$ is a linear transformation, we can take the ordered basis $\mathcal{B} := \{e_{1},e_{2}\}$ and the direction $v = (v_{1},v_{2})$ in order to conclude that
\begin{align*}
D_{v}f(p) = f'(p)v & = f'(p)(v_{1}e_{1} + v_{2}e_{2})\\\\
& = v_{1}f'(p)e_{1} + v_{2}f'(p)e_{2}\\\\
& = v_{1}D_{e_{1}}f(p) + v_{2}D_{e_{2}}f(p)\\\\
& = v_{1}\frac{\partial f}{\partial x}(p) + v_{2}\frac{\partial f}{\partial y}(p)
\end{align*}
Solution
Now it remains to apply the mentioned result to the particular case where $f(x,y) = x^{2} - 3xy + 4y^{2}$.
Can you take it from here?
