The wave equation 1D with constant density is defined as:
\begin{equation} \frac{\partial^2 U}{\partial t ^2} = V^2 \frac{\partial^2 U}{\partial x ^2} \label{eqa} \end{equation}
And the implicit difference schema:
\begin{equation} U_j^{n+1} = \left( \frac{\Delta t V_j^n}{\Delta s} \right) ^2 \left( U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n - U_j^{n-1} \label{eq:b} \end{equation}
Where $\Delta x = \Delta s $.
Von Newman stability analysis.
Since $ u(x,t) = e^{i(wt+kx)} $ is a solution for the discrete schema above can be written as:
\begin{eqnarray} u(x_j,t_n) &=& e^{i(wt_n+kx_j)} \nonumber \\ &=& e^{iwn\Delta t} e^{ikj\Delta s} \nonumber \\ &=& \epsilon^n e^{ikj\Delta s} \label{eq:c} \end{eqnarray}
Where $\epsilon = e^{iw \Delta t} $ and $ i = \sqrt{-1} $ . To maintain stability we should make sure that $ \epsilon \leq 1 $ not growing exponentially with increasing time steps. So applying the previous equation in the implicit differences schema we can analyse the growth.
\begin{multline} U_j^{n+1} = \left( \frac{\Delta t V_j^n}{\Delta s} \right) ^2 \left( U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n - U_j^{n-1} \\ \epsilon^{n+1} e^{ikj\Delta s} = r^2 \left( \epsilon^{n+1} e^{ik(j+1)\Delta s} - 2 \epsilon^{n+1} e^{ikj\Delta s} + \epsilon^{n+1} e^{ik(j-1)\Delta s} \right) \\ + 2 \epsilon^n e^{ikj\Delta s} - \epsilon^{n-1} e^{ikj\Delta s} \\ 1 = r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 2 \epsilon^{-1} - \epsilon^{-2} \\ \epsilon^{-2} - 2 \epsilon^{-1} - r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 1 = 0 \\ \epsilon^{-2} - 2 \epsilon^{-1} - 4 r^2 \sin^2 \left( k \Delta s / 2 \right) + 1 = 0 \label{eq:d} \end{multline}
Where $ r = \left( \frac{\Delta t V_j^n}{\Delta s} \right) $
At the previous equation we can substitute $ \phi = \epsilon^{-1} $ turning it to a second degree bellow:
\begin{eqnarray} \phi^2 - 2 \phi -4 r^2 \sin^2 \left(k\Delta s/2 \right) +1 = 0 \nonumber \\ \phi^2 - 2 \phi + c = 0 \nonumber \end{eqnarray}
With: $$ c = 1 -4 r^2 \sin^2 \left(k\Delta s/2 \right) $$
Has roots $ \phi^{'} $ and $ \phi^{''} $ as bellow. Replacing $c$ also.
\begin{eqnarray} &=& \frac{2 \pm \sqrt{4 - 4c}}{2} = 1 \pm \sqrt{1-c} \nonumber \\ &=& 1 \pm \sqrt{4 r^2 \sin^2 \left(k\Delta s/2 \right) } \nonumber \\ &=& 1 \pm 2 r \sin \left(k\Delta s/2 \right) \nonumber \\ \phi^{'} &=& 1 + 2 r \sin \left(k\Delta s/2 \right) \\ \phi^{''} &=& 1 - 2 r \sin \left(k\Delta s/2 \right) \end{eqnarray}
Going back to $\epsilon$ we get :
\begin{eqnarray} \epsilon^{'} &=& \frac{1}{1 + 2 r \sin \left(k\Delta s/2 \right)} \\ \epsilon^{''} &=& \frac{1}{1 - 2 r \sin \left(k\Delta s/2 \right) } \end{eqnarray}
But with $r \sin \left(k\Delta s/2 \right) \to 0$ we get that $\epsilon^{''} \to \infty $
Should it not be unconditionally stable? What is wrong in here?
Update: After awesome answer by Jitse, fixing the signal we get:
\begin{eqnarray} \phi^{'} &=& 1 + i 2 r \sin \left(k\Delta s/2 \right) \\ \phi^{''} &=& 1 - i 2 r \sin \left(k\Delta s/2 \right) \end{eqnarray}
Analyzing the modulus, since $ \phi^{'} $ and $ \phi^{''} $ are conjugate pairs of the same complex number the modulus is the same for both. We get:
$$ \| \phi^{'} \| = \| \phi^{''} \| = \sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}$$
So since $ \| \phi \| \geq 1 \ \ \forall \ r, \ k,\ \Delta s $ than
$$ \| \epsilon \| = \frac{1}{ \sqrt{1 + 4 r ^2 \sin ^2 \left(k\Delta s/2 \right)}} $$
$$ \| \epsilon \| \leq 1 \ \ \forall \ r, \ k,\ \Delta s$$
Awesome, unconditionally stable!!
Your method for Neumann stability analysis is correct, and indeed I would expect the scheme to be unconditionally stable.
I do not understand the step where you introduce the sine. I think that $$ e^{ix} - 2 + e^{-ix} = -4\sin^2(\tfrac12x). $$ The minus sign seems to have been lost in your computation.