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.
Actually, Schroedinger equation
$$ -i\hbar \partial_t \psi+ \left[ -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)^2+qV \right]\psi=0\tag{0}$$
under the gauge transformations
$$ \psi \rightarrow \psi'= e^{iq\Lambda/\hbar} \psi$$
$$ \vec{A} \rightarrow \vec{A}' = \vec{A} + \nabla \Lambda$$
$$ V \rightarrow V'= V - \partial_t\Lambda$$
does not remain invariant, but the left-hand side of (0) gives rise to
$$ -i\hbar \partial_t \psi'+ \left[ -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}'\right)^2+qV' \right]\psi'=
e^{iq\Lambda/\hbar}\left\{-i\hbar \partial_t \psi+ \left[ -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)^2+qV \right]\psi\right\}\:.$$
In summary, since $e^{iq\Lambda/\hbar}\neq 0$,
$\qquad\quad$ gauge transformed quantities satisfy Schroedinger equation if untranformed quantities do.
To prove it, avoid brute force computations as yours which give rise to unavoidable mistakes almost certainly and go on as follows. First rewrite the initial equation as
$$ -\left[i\hbar \partial_t -qV \right]\psi -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)^2\psi=0\tag{1}$$
Next notice that, under the transformations, we have
$$\left[i\hbar \partial_t -qV' \right]\psi' =
\left[i\hbar \partial_t -q(V - \partial_t\Lambda)\right]e^{iq\Lambda/\hbar}\psi =
e^{iq\Lambda/\hbar}\left[i\hbar \partial_t -qV \right]\psi $$
and
$$\left(\frac{\hbar}{i}\nabla-q\vec{A}'\right)\psi' = \left(\frac{\hbar}{i}\nabla-q(\vec{A} +\nabla \Lambda)\right)e^{iq\Lambda/\hbar}\psi = e^{iq\Lambda/\hbar}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)\psi$$
so that, iterating the second result
$$\left(\frac{\hbar}{i}\nabla-q\vec{A}'\right)^2\psi' = e^{iq\Lambda/\hbar}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)^2\psi\:.$$
Putting all together, under the action of gauge transformations, (1) becomes
$$ -\left[i\hbar \partial_t -qV' \right]\psi' -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}'\right)^2\psi' =
e^{iq\Lambda/\hbar}\left\{-\left[i\hbar \partial_t -qV \right]\psi -\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-q\vec{A}\right)^2\psi\right\}=0$$
as wanted.
Best Answer
If $\Lambda$ is a function of the position operator, you can't just treat $\Lambda$ as a scalar function. In fact, what you are obtaining is $\tilde{\mathbf{P}} = \mathbf{P} - i q [\mathbf{P},\Lambda(\mathbf{R})]$.
For analytic functions you can always expand $\Lambda$ as a power series in $\mathbf{R}$ and compute the commutator with $\mathbf{P}$ for each power of $\mathbf{R}$ using the canonical commutation relations and the product rule for operators. You will obtain $[\mathbf{P},\Lambda(\mathbf{R})] = - i \Lambda'(\mathbf{R})$. Using this you get the result you were looking for.