Why my 4-divergence term added to a Lagrangian modifies the equation of motion? I take this Lagrangian: 
$$\mathcal{L}=\mathcal{L}_0+\partial_\alpha f(\phi, \partial_\mu \phi).$$
In this topic Does a four-divergence extra term in a Lagrangian density matter to the field equations? , it is said that any 4-divergence term added to a Lagrangian doesn't modifies the equation of motion.
In my example I add $\partial_\alpha f(\phi, \partial_\mu \phi)$ to $\mathcal{L}_0$ (it is not a 4-divergence but the mechanics behind is exactly the same). And I remark that it can modify the equation of motion if $f$ contains time derivatives of $\phi$. So I don't understand.
I Write the infinitesimal variation of action to $\mathcal{L}$:
$$ \delta S = \int d^4x ~ \delta \mathcal{L}, $$
$$ \delta S = \int d^4x ~ [ \frac{\partial \mathcal{L}_0}{\partial \phi} \delta \phi +  \frac{\partial \mathcal{L}_0}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi) + \partial_\alpha [\frac{\partial f}{\partial \phi} \delta \phi + \frac{\partial f}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi)] ~ ].$$
As usual, I know that : $\delta(\partial_\mu \phi)=\partial_\mu \delta(\phi)$. Thus I can integrate by parts: 
$$ \delta S = \int d^4x ~ [ \frac{\partial \mathcal{L}_0}{\partial \phi}  - \partial_\mu  \frac{\partial \mathcal{L}_0}{\partial (\partial_\mu \phi)} )\delta \phi + \int d^4x ~ \partial_\mu[\frac{\partial \mathcal{L}_0}{\partial (\partial_\mu \phi)} \delta \phi] + \int d^4x ~ \partial_\alpha [\frac{\partial f}{\partial \phi} \delta \phi + \frac{\partial f}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi)].$$
We have: 
$$ \int d^4x ~ \partial_\mu[\frac{\partial \mathcal{L}_0}{\partial (\partial_\mu \phi)} \delta \phi] = \int d^3x ~ [\frac{\partial \mathcal{L}_0}{\partial (\partial_\mu \phi)} \delta \phi]_{x_i^{-}}^{x_i^{+}}=0.$$
Indeed, $\delta \phi=0$ on the boundaries by hypothesis ($x_i^{+}=+\infty$ for spatial coordinates and $t_f$ for time).
We also have:
$$ \int d^4x ~ \partial_\alpha [\frac{\partial f}{\partial \phi} \delta \phi + \frac{\partial f}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi)]= \int d^3x ~ [\frac{\partial f}{\partial \phi} \delta \phi + \frac{\partial f}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi)]_{x_i^{-}}^{x_i^{+}}=\int d^3x ~ [\frac{\partial f}{\partial (\partial_\mu \phi)} \delta(\partial_\mu \phi)]_{x_i^{-}}^{x_i^{+}}.$$
** And here is my problem **.
The fact $\delta \phi(x_i^{+})=\delta \phi(x_i^{-})=0$ doesn't implicate that $\partial_\mu \delta \phi(x_i^{+})=\partial_\mu \delta \phi(x_i^{-})=0$.
To be more precise, it could be true if $x_i^{+}=-x_i^{-}=+\infty$(*) but if I take the time coordinates, I have $x_i^{+}=t_f$. So it is at least not true for $\mu=t$.
Thus the extra term $\partial_\alpha f(\phi, \partial_\mu \phi)$ modifies the extremality of the action. Thus I will not have the same equation of motion.
But in this topic : Does a four-divergence extra term in a Lagrangian density matter to the field equations? the book of the author says that any four divergence doesn't affect the equation of motions.
But we've seen here (if I made no mistake which is not sure at all) that if the extra term is a total derivative that contains time derivatives of the field it can change the equations of motion.
Where am I wrong?

(*) : it is true because we ask $\phi$ to go to zero at infinity, so we only allow variations of $\phi$ that vanish at infinity (else we would end up with $\phi+\delta \phi$ not integrable). And as $(x,y,z) \mapsto \delta \phi(x,y,z,t)$ goes to $0$ at infinity, all its derivative also.
 A: The correct statement is that a boundary term (BT) in the action (or equivalently, a total divergence term in the Lagrangian density) does not change the functional/variational derivative if both the old and the new functional derivatives exist. Pay attention to the important word if in the previous sentence: This does not exclude the possibility that a functional/variational does not exist. 
In order for functional derivatives to exist, it is necessary to impose adequate boundary conditions (BCs). A boundary/total divergence term may change the adequate set of BCs.
In OP's example, he has correctly observed that Dirichlet BCs are not enough to remove BTs in the variation.
To summarize: OP has not shown that 2 different sets of Euler-Lagrange equations exist, cf. the title question (v6). Only that some choices of BTs & BCs may make the variational problem ill-defined.
For the point mechanical case, see also this Phys.SE post. The field theoretical case is a straightforward generalization.
A: On the boundary, $\delta \phi(x) = 0 \implies \delta \left(\partial_{\mu} \phi(x)\right) = 0$. 
Think in terms of the one-dimensional variational principle. In this case, one finds the equivalence $\delta \phi(x) = \delta \left(\dot{\phi}(x) dt\right) = \delta \left(\dot{\phi}(x)\right) dt$. Thus when one takes $\delta \phi(x) = 0$ on the boundary, we obtain immediately $\delta \left(\dot{\phi}(x)\right) = 0$ as well.
This holds for any variational principle with the given boundary condition in any dimension. I hope this solves your confusion.
