How to calculate functional derivative correctly?

Once nice way to calculate functional derivatives is to use the concept of the Gateaux derivative as follows:

$$\frac{d}{d\epsilon}S[\phi+\epsilon \eta]\bigg|_{\epsilon=0} = \int d^4x\frac{\delta S}{\delta \phi} \eta$$

In your case, $$S[\phi+\epsilon \eta]= \int d^4x \ \bigg\{\frac{1}{2}\big((\partial \phi)^2 + 2\epsilon (\partial_\mu\phi)(\partial^\mu\eta) + \epsilon^2(\partial \eta)^2\big)$$ $$ + \frac{1}{2}m^2\big(\phi^2+2\epsilon \eta\phi+\epsilon^2\eta^2\big)+ V(\phi+\epsilon \eta)- J(\phi+\epsilon\eta)\bigg\}$$ Differentiating and setting $\epsilon$ to zero yields $$ \frac{d}{d\epsilon}S[\phi+\epsilon\eta]\bigg|_{\epsilon=0} = \int d^4x \bigg\{(\partial_\mu\phi)(\partial^\mu\eta) + m^2\phi \eta + V'(\phi)\eta - J\eta\bigg\}$$ We can cast this into the desired form by integrating by parts, yielding $$\frac{d}{d\epsilon}S[\phi+\epsilon\eta]\bigg|_{\epsilon=0} = \int d^4x\bigg\{-\partial^2\phi + m^2\phi + V'(\phi) - J\bigg\} \eta$$ We can therefore read off $$\frac{\delta S}{\delta \phi} = -\partial^2\phi + m^2 \phi+ V'(\phi) - J$$

Here is a second way to see the correct result for taking the functional derivative of the spacetime derivative of the field, which I hope will be helpful.

Recall that the definition of the functional derivative is $$ \frac{\delta\phi(y)}{\delta\phi(x)}=\delta(y-x) .$$ You know that Dirac deltas are distributions. That is, you should always think of them living under an integral with some test function. So the above definition should really be thought of as $$ \frac{\delta}{\delta\phi(x)} \int \phi(y) f(y)\, dy = \int \delta(y-x) f(y)\, dy=f(x)$$ for some arbitrary function $f(y)$.

Now suppose instead you have the spacetime derivative of $\phi$. $$ \frac{\delta}{\delta\phi(x)} \int \partial\phi(y) f(y)\, dy $$ To understand what this means, just integrate by parts. $$ -\frac{\delta}{\delta\phi(x)} \int \phi(y)\, \partial f(y)\, dy= -\int \delta(y-x)\, \partial f(y)\, dy =-\partial f(x) $$ But this is exactly the definition of how the derivative of the Dirac delta is supposed to act on an arbitrary test function. Informally, you can just integrate by parts back to get $$ -\int \delta(y-x)\, \partial f(y)\, dy = \int \partial\delta(y-x)\, f(y)\, dy .$$ Pulling this result out of its nice safe integral home, we can write the definition $$\frac{\delta}{\delta\phi(x)}\partial\phi(y) = \delta'(y-x).$$

Applying this definition to your problem gives the desired result. In shorthand, we can say that the functional derivative just "moves past" the spacetime derivative on the field ($\delta\partial\phi=\partial\delta\phi$), so that it acts "as you'd expect" on the two factors of $\partial\phi$ and gives you the factor of 2 that you need.

The safest way to compute the functional derivative is to use the following prescription:

\begin{equation} S[\phi + \delta \phi] = S[\phi] + \int {\rm d}^4 x \frac{\delta S}{\delta \phi}\delta \phi + O(\delta \phi^2) \end{equation}

In other words, add a small perturbation to the field, and manipulate the action so it has the form of an integral times the variation (ignoring terms higher than linear order in the variation). Then the part of the integrand multiplying the variation is the functional derivative.

Here's how to apply this in your example.

We start with the action (I am going to absorb the mass term into the potential since it doesn't really make any difference for this calculation) \begin{equation} S[\phi] = \int {\rm d}^4 x \left(\frac{1}{2}(\partial \phi)^2 + V(\phi) + \phi J \right) \end{equation}

Then we add a perturbation to the field and only keep terms to first order

\begin{equation} S[\phi+\delta \phi] = S[\phi] + \int {\rm d}^4 x \left(\partial_\mu \phi \partial^\mu \delta \phi + \frac{\partial V}{\partial \phi}\delta \phi + \delta \phi J \right) + O(\delta \phi^2) \end{equation}

Then we do an integration by parts on the kinetic term so that we remove the derivative from the variation. This leads to

\begin{equation} S[\phi+\delta \phi] = S[\phi] + \int {\rm d}^4 x \left[\left(-\square \phi + \frac{\partial V}{\partial \phi}+ J \right) \delta \phi \right] + O(\delta \phi^2) \end{equation}

Comparing with the definition above, we see that

\begin{equation} \frac{\delta S}{\delta \phi} = -\square \phi + \frac{\partial V}{\partial \phi} + J \end{equation}