Hmmm... First of all, if phi is a function of q and epsilon (i.e. \phi=\phi(q,\epsilon) ), I think in general phi-dot should be a function of q, epsilon and their time derivatives (i.e. \dot{\phi}=\dot{\phi}(q,\dot{q},\epsilon,\dot{\epsilon})
Secondly, you need to be more mindful about which derivatives are partials and which are normal derivs...for example, since you are taking the normal derivative of your Lagrangian you should have:
\frac{dL}{d\epsilon} = \frac{\partial L}{\partial \phi}\frac{d \phi}{d \epsilon} + \frac{\partial L}{\partial \dot{\phi}}\frac{d \dot{\phi}}{d \epsilon}
The way you had it written, it would simplify under the usual rules for partial derivatives to \frac{dL}{d\epsilon} = \frac{\partial L}{\partial \epsilon}+ \frac{\partial L}{\partial \epsilon}=2\frac{\partial L}{\partial \epsilon} which is clearly nonsense.
Thirdly, if \dot{\phi} = \frac{d\phi}{dt}= \frac{d\phi}{dq}\frac{dq}{dt} = \frac{d\phi}{dq}\dot{q}, then shouldn't the product rule actually give you: \frac{d \dot{\phi}}{d \epsilon} = \frac{ d^2 \phi}{d q d \epsilon}\dot{q}+ \frac{d\phi}{dq} \frac{d\dot{q}}{d \epsilon}
Finally, I think the best way to go about things is to use the fact that \phi=\phi(q,\epsilon) and \dot{\phi}=\phi(q,\dot{q},\epsilon,\dot{\epsilon}) to compute the derivatives \frac{d \phi}{d \epsilon} and \frac{d \dot{\phi}}{d \epsilon} in terms of \frac{d q}{d \epsilon} and \frac{d \dot{q}}{d \epsilon} using the chain rule in the same way you did to obtain \frac{dL}{d\epsilon} = \frac{\partial L}{\partial \phi}\frac{d \phi}{d \epsilon} + \frac{\partial L}{\partial \dot{\phi}}\frac{d \dot{\phi}}{d \epsilon}.
You should find that you are able to cancel your \partial \phi and \partial \dot{\phi} derivs