Calculus method
The Taylor series of a function of $d$ variables is as follows:
$$f(\mathbf{x} + \mathbf{y}) = \sum_{n_1=0}^{\infty}\ldots\sum_{n_d=0}^{\infty}
\frac{(y_1-x_1)^{n_1} \ldots (y_d - x_d)^{n_d}}{n_1!\ldots n_d!}
\left. \left(
\partial_1^{n_1}\ldots \partial_d^{n_d} \right) f \right|_{\mathbf{x}}.$$
where $\partial_i^{n_i} f$ means "the $n_i^{\text{th}}$ order partial derivative of $f$ with respect to the $i^{\text{th}}$ coordinate".
Consider only the terms where $\sum_{i=1}^{\infty} n_i = 1$.
These are the terms for which there is only one derivative being taken, a.k.a. the "linear terms".
Keeping the $0^{\text{th}}$ order term and the linear terms gives
$$f(\mathbf{x} + \mathbf{y}) = f(\mathbf{x})+\sum_{i=1}^{d}(y_i - x_i) \left. (\partial_i f \right)|_{\mathbf{x}} . \qquad (*)$$
Define the displacement $\epsilon = \mathbf{y} - \mathbf{x}$. Then
$$\epsilon \cdot \nabla = \sum_{i=1}^d (y_i - x_i) \partial_i$$
by definition of what $\cdot$ means. Therefore, keeping only up to linear terms, $(*)$ becomes
$$f(\mathbf{y} - \mathbf{x}) = f(\mathbf{x}) + \left. (\epsilon \cdot \nabla)\right|_{\mathbf{x}}f . $$
So you see, your formula is correct, just in a particular notation you may not be used to.
Algebraic method
For the sake of typing less, I'm going to explain this in one dimension, but nothing in the actual calculation is restricted to one dimension.
The state vector can be thought of as a linear combination of eigenvectors of the position operator
$$|\Psi\rangle = \int_x \Psi(x) |x\rangle \qquad (1)$$
where $\hat{x}|x\rangle = x|x\rangle$.
Bringing in a bra $\langle y |$ from the left and using $\langle y | x \rangle = \delta(x-y)$ in (1) gives
$$\langle y | \Psi \rangle = \Psi(y) . \qquad (2)$$
Combining (1) and (2) gives
$$|\Psi\rangle = \int_x |x\rangle \langle x| \Psi \rangle \qquad (3)$$
which shows that $\int_x |x\rangle \langle x | = \text{Identity}$.
Let's really understand what all this means:
Eq. (1) just says that you can write a vector as a linear combination of basis vectors.
In this case, $\Psi(x)$ are the coefficients in the linear combination.
In other words, the wave function is just the coefficients of a linear combination expansion written in the position basis.
Eq. (2) makes this explicit by showing that the inner product of a position basis vector $|y\rangle$ with $|\Psi\rangle$ is precisely $\Psi(y)$.
Eq. (3) just expresses a neat way of expressing the identity operator in a way which we will find very useful.
Note that this works with any basis.
For example,
$$\int_p |p\rangle \langle p | = \text{identity} .$$
We can use this to express a position eigenvector in terms of momentum eigenvectors,
$$|x\rangle = \int_p |p\rangle \langle p | x \rangle = \int_p e^{-i x p/\hbar}|p\rangle ,$$
where we've used the fact that $\langle x | p \rangle = e^{i p x/\hbar}$ [1].
Now let's get back to the question at hand. Define $|\Psi'\rangle \equiv e^{i \epsilon \hat{p}/\hbar}|\Psi\rangle$.
Let us evaluate $\Psi'(y) \equiv \langle y|\Psi'\rangle$:
\begin{align}
\Psi'(y) &= \langle y | \Psi' \rangle \\
&= \langle y| e^{i \epsilon \hat{p}/\hbar} |\Psi \rangle \\
&= \langle y| e^{i \epsilon \hat{p}/\hbar} \int_x \Psi(x) |x\rangle \\
&= \int_x \Psi(x) \langle y | e^{i \epsilon \hat{p}/\hbar} |x \rangle .
\end{align}
To proceed we need to compute $e^{i \epsilon \hat{p}/\hbar} |x \rangle$.
We can do this using our expression for the identity in terms of momentum states,
$$
\begin{align}
e^{i \epsilon \hat{p}/\hbar} |x \rangle
&= e^{i \epsilon \hat{p}/\hbar} \left( \int_p |p\rangle \langle p | \right) | x \rangle \\
&= \int_p e^{i \epsilon \hat{p} / \hbar} |p\rangle \langle p | x \rangle \\
&= \int_p e^{i \epsilon p / \hbar} |p\rangle \langle p | x \rangle \\
&= \int_p e^{i \epsilon p /\hbar} |p\rangle e^{-i p x/\hbar}\\
&= \int_p e^{-i(x - \epsilon) p / \hbar} |p\rangle \\
&= |x - \epsilon \rangle
\end{align}$$
This is actually the thing you should remember about this question:
$$e^{i\epsilon\hat{p}/\hbar}|x\rangle = |x - \epsilon \rangle .$$
This is an excellent equation because it tells you something about how an the $e^{i \epsilon \hat{p} /\hbar}$ operator changes a position eigenvector without having to write it out in a particular basis.
Now that we have that, we can go back to what we were trying to compute,
$$
\begin{align}
\Psi'(x) &= \int_x \Psi(x) \langle y | e^{i \epsilon \hat{p}/\hbar} |x \rangle \\
&= \int_x \Psi(x) \langle y | x - \epsilon \rangle \\
&= \int_x \Psi(x) \delta(y - (x - \epsilon)) \\
&= \Psi(y + \epsilon) .
\end{align}
$$
This is the equation you wanted to prove.
If I messed up any minus signs I hope someone will edit :)
[1] I'm not proving that, and if you want to know why that's true it would make a good question on this site.
