If we use the linear Lagrange basis for finite element method the semi-discrete solution is $$ u_h(x,t) = \sum_{i=1}^n u_i(t)\phi_i(x).$$
For the homogenous Neumann boundary, the solution $u_h$ have to satisfy the following condition $$\frac{\partial u_h}{\partial x} = \frac{u_1 - u_2}{h} =0,$$ which is invalid for some cases. In some papers, they said that they use a linear piecewise polynomial. How they apply the basis with this problem? I just doubt it.
To illustrate the idea, let's consider the Helmholtz equation $-\Delta u + u = f$ with homogeneous Dirichlet BCs compared to Homogenous Neumann BCs. In both cases the linear and bilinear form are the same, and the weak problem is to find $u\in V$ such that
$$\int \nabla u\cdot\nabla v + uv = \int fv\quad\forall v\in V,$$
where $V=H^1_0$ in the Dirichlet case, or $V=H^1$ in the Neumann case. Since the bilinear and linear forms are the same, we can see that the BC enforcement comes from the choice of space. In the Dirichlet case the BCs are imposed on the space directly, whereas for the Neumann case the BC is imposed weakly in the bilinear form (a function u had to satisfy the Neumann BC for the weak formulation to be valid when the test functions do not vanish on the boundary).
All of this carries over to the finite element setting analogously, and the difference is in the Dirichlet case one would use the Lagrange space with functions that vanish on the boundary, whereas in the Neumann case the finite element space is larger, and has functions that are nonzero on the boundary.