One often comes across stability regions when looking at explicit and implicit Euler's method (for $\dot{x}=\lambda x$). But I have never come across such region for Verlet, say for the ODE $\ddot{x}=\lambda x$: $$ x_{n+1} = 2x_n - x_{n-1} + h^2\lambda x_n $$ So I thought I would try this my self. After initially writing out quite a few terms to find a pattern I found the following expression. First I will define the "growth" parameter $\alpha$ as $$ \alpha = 2 + h^2 \lambda $$ $$ x_{k+1} = x_1 \sum^{\lfloor \frac{k}{2}\rfloor}_{m=0} \left(\frac{(-1)^m \alpha^{k-2m}}{m!} \prod^{2m-1}_{i=m}(k-i)\right) - x_0 \sum^{\lfloor \frac{k-1}{2}\rfloor}_{m=0} \left(\frac{(-1)^m \alpha^{k-1-2m}}{m!} \prod^{2m-1}_{i=m}(k-1-i)\right) $$ For instance for $k=10$ this yields: $$ x_{11} = x_1 \left(\alpha^{10} - 9\alpha^8 + 28\alpha^6 - 35\alpha^4 + 15\alpha^2 - 1\right) - x_0 \left(\alpha^9 - 8\alpha^7 + 21\alpha^5 - 20\alpha^3 + 5\alpha\right) $$ However I am wondering if these terms would converge or diverge depending on the value of $\alpha$. I have no idea how to approach this since there is product within summation and I already found it quite an achievement when I found this expression.
This might be a related question, since it also looks a the stability behavior of Verlet, but I would like to known whether these expressions would converge or diverge for large $k$ depending on $\alpha$.
Edit: The answer of Did did convince me that approaching this problem with matrices would be the easier way to go, but Did did not proof that what values for $c_1$, $c_2$, $c_3$ and $c_4$ are, so it might be possible that the limit of $n$ to infinity would not blow up.
First I will, for clarity sake, define the matrix notation of the problem which has to solved,
$$ \begin{bmatrix}x_{n+1} \\ x_{n}\end{bmatrix} = A \begin{bmatrix}x_{n} \\ x_{n-1}\end{bmatrix} $$
$$ A = \begin{bmatrix}\alpha & -1 \\ 1 & 0\end{bmatrix}. $$
Thus by multiplying the initial positions vector $n$ times by $A$ will allow you to find the positions vector $n$ steps further,
$$ \begin{bmatrix}x_{n+1} \\ x_{n}\end{bmatrix} = A^n \begin{bmatrix}x_{1} \\ x_{0}\end{bmatrix}. $$
The values for those constants can be derived with the help of eigen decomposition of a matrix $A$,
$$ A = PDP^{-1}, $$
Matrix $D$ is a diagonal matrix with each element of the diagonal containing one of the eigenvalues of $A$, which are
$$ \lambda_1 = \frac{\alpha + \sqrt{\alpha^2 - 4}}{2}, $$
$$ \lambda_2 = \frac{\alpha - \sqrt{\alpha^2 - 4}}{2}. $$
By expressing $\lambda_1$ in $\lambda_2$ it can be shown that
$$ \lambda_1 = \lambda_2^{-1}. $$
In order to keep the rest of the equations a little less bulky I will use $\lambda$ for the first eigenvalue and thus the second eigenvalues will be $\lambda^{-1}$, instead of the parameter $\alpha$. The resulting matrices $D$, $P$ and $P^{-1}$ are therefore
$$ D = \begin{bmatrix}\lambda & 0 \\ 0 & \lambda^{-1}\end{bmatrix}, $$
$$ P = \begin{bmatrix}\lambda & \lambda^{-1} \\ 1 & 1\end{bmatrix}, $$
$$ P^{-1} = \frac{1}{\lambda^2 - 1}\begin{bmatrix}\lambda & -1 \\ -\lambda & \lambda^2\end{bmatrix}. $$
By substituting these matrices for $A$ in the expressing for the positions vector $n$ steps further yields
$$ \begin{bmatrix}x_{n+1} \\ x_{n}\end{bmatrix} = P D^n P^{-1} \begin{bmatrix}x_{1} \\ x_{0}\end{bmatrix} = \frac{1}{\lambda^2 - 1} \begin{bmatrix}\lambda & \lambda^{-1} \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda^n & 0 \\ 0 & \lambda^{-n}\end{bmatrix} \begin{bmatrix}\lambda & -1 \\ -\lambda & \lambda^2\end{bmatrix} \begin{bmatrix}x_{1} \\ x_{0}\end{bmatrix}. $$
When multiplying this all out you will get
$$ x_{n+1} = x_1 \frac{\lambda^{n+2} - \lambda^{-n}}{\lambda^{2} - 1} + x_0 \frac{\lambda^{1-n} - \lambda^{n+1}}{\lambda^{2} - 1}. $$
I am not sure how to proceed here, since $\lambda$ can be complex, when $2<\alpha<-2$, in which case its magnitude will also always be one. I also did some numerical calculations with MATLAB and it seemed to be stable for $2<\alpha<-2$, however I do not know how to proof this. And stability regions also include complex values for $h$.
This is really the same question since the recursion is $v_{n+1}=Av_n$ where $v_n$ is the column vector $v_n=\begin{pmatrix}x_{n+1}\\ x_{n}\end{pmatrix}$ and $A$ is the $2\times2$ matrix $A=\begin{pmatrix}a&-1\\1&0\end{pmatrix}$. Thus $v_n=A^nv_0$ for every $n$ and, if the eigenvalues $\lambda_1$ and $\lambda_2$ of $A$ are distinct, there exists some constants $(c_1,c_2,c_3,c_4)$ such that $x_n=(c_1\lambda_1^n+c_2\lambda_2^n)x_0+(c_3\lambda_1^n+c_4\lambda_2^n)x_1$ for every $n\geqslant0$.
If $|\lambda_1|\lt1$ and $|\lambda_2|\lt1$, then, for every initial condition $(x_1,x_0)$, $x_n\to0$ when $n\to\infty$. Finally, note that $\lambda_1$ and $\lambda_2$ are the roots of the characteristic polynomial $\chi_A(x)=\det(xI-A)=x^2-ax+1$ hence $\lambda_1\lambda_2=1$ and it is impossible that $|\lambda_1|\lt1$ and $|\lambda_2|\lt1$.
Note that $a=2+h^2\lambda$ hence $\chi_A(1)=\chi'_A(1)=2-a=-h^2\lambda$. If $h^2\lambda\gt0$, then $\chi_A(1)=\chi'_A(1)\lt0$ hence one eigenvalue is real and larger than $1$ and there is explosion.