Finding the general term of two related recurrence relations I'm trying to find the general term of the recurrence relations

$\quad a_{n+1}=a_n+\text hb_n$
$\quad b_{n+1}=b_n-\text ha_n $
$\quad a_0=0, \quad b_0=1$

I tried finding the terms, $a_1, \space a_2, \cdots$ so I can find a formula for the general term but it was too long and and I couldn't even find the solution. Is there a better way I can approach it?

To be specific, I found these relations while using the explicit Euler method for solving ODE's on a second order DE. I'm will be working on this a lot and I was wondering if there is a better way to find the general term for these recurrence relations other than checking term by term until I can find a pattern?

 A: Assuming $h\neq0$ is constant, from the first equation we get: $b_n=\frac1h\left(a_{n+1}-a_n\right)$. Substituting into the second, we have:
$$b_{n+1}=\frac1h\left(a_{n+2}-a_{n+1}\right)=\frac1h\left(a_{n+1}-a_n\right)-ha_n \Rightarrow\\0= a_{n+2}-a_{n+1}-a_{n+1}+a_n+h^2a_n=a_{n+2}-2a_{n+1}+(h^2+1)a_n\Rightarrow\\
a_{n}-2a_{n-1}+(h^2+1)a_{n-2}=0, \hspace{10pt} a_0=0, a_1=0+h=h$$
Now, following my explanation here, we will solve the last recurrence relation:
The polynomial is $t^2-2t+h^2+1$, hence $t_{1,2}=1\pm ih$.
So $a_n=\alpha(1+ih)^n+\beta(1-ih)^n$. Using the values of $a_0,a_1$, we can find that $\alpha=-\frac i2$ and $\beta=\frac i2$.
So $$\begin{array}{l}
a_n=\frac i2\left((1-ih)^n-(1+ih)^n\right)\\
\begin{align*}b_n=&\frac1h\left(a_{n+1}-a_n\right)=\frac i{2h}\left((1-ih)^{n+1}-(1+ih)^{n+1}-(1-ih)^n+(1+ih)^n\right)\\
=&\frac i{2h}\left(-ih(1-ih)^n-ih(1+ih)^n\right)=\frac12\left((1-ih)^n+(1+ih)^n\right)\end{align*}
\end{array}$$
A: You can try working as follows, which gives a longer recursion in one variable.$$a_{n+2}=a_{n+1}+hb_{n+1}=a_{n+1}+hb_n-h^2a_n$$$$=a_{n+1}+(a_{n+1}-a_n)-h^2a_n$$
And on tidying up this becomes $$a_{n+2}=2a_{n+1}-(1+h^2)a_n$$
A: Introducing the line vectors $x_n=(a_n,b_n)$, one sees that $x_0=(0,1)$ and $x_{n+1}=x_nM$ with 
$$
M=\begin{pmatrix}1 & -h\\ h & 1\end{pmatrix}.
$$ 
The characteristic polynomial of $M$ is 
$$
\chi_M(\lambda)=\det(\lambda I-M)=\lambda^2-2\lambda+1+h^2,
$$ 
whose roots are $\lambda_\pm=1\pm\mathrm ih=r\mathrm e^{\pm\mathrm i\theta}$ with $r=\sqrt{1+h^2}$, $r\cos\theta=1$, $r\sin\theta=h$, and a left-eigenvector for $\lambda_\pm$ is $v_\pm=(1,\pm\mathrm ih)$, that is, $v_\pm M=\lambda_\pm v_\pm$.
If $x_0=t_+v_++t_-v_-$ for some $t_\pm$, then $x_n=t_+\lambda_+^nv_++t_-\lambda_-^nv_-$. Solving for $t_\pm$ yields $0=t_++t_-$, $1=t_+\mathrm ih-t_-\mathrm ih$, hence $t_+=-t_-=\frac1{2\mathrm ih}$. 
Thus, $a_n=x_n\cdot(1,0)^T=t_+\lambda_+^n+t_-\lambda_-^n=\frac1{2\mathrm ih}r^n(\mathrm e^{n\mathrm i\theta}-\mathrm e^{-n\mathrm i\theta})$, that is,
$$
a_n=\frac1hr^n\sin(n\theta),\quad r=\sqrt{1+h^2},\quad \sin(\theta)=\frac{h}r,\quad |\theta|\leqslant\frac\pi2.
$$
For every $n\geqslant1$, an equivalent formulation is
$$
a_n=r^{n-1}U_{n-1}(1/r),\qquad r=\sqrt{1+h^2},
$$
where $U_n$ is the $n$th Chebyshev polynomial of the second kind. Likewise, for every $n\geqslant0$,
$$
b_n=r^n\cos(n\theta)=r^nT_n(1/r),\qquad r=\sqrt{1+h^2},
$$
where $T_n$ is the $n$th Chebyshev polynomial of the first kind.
A: Set up generating functions $A(z) = \sum_{n \ge 0} a_n z^n$ and similarly $B(z)$. Multiply the recurrences by $z^n$, sum over $n \ge 0$, and by recognizing some sums you get:
\begin{align}
\frac{A(z) - a_0}{z} &= A(x) + h B(z) \\
\frac{B(z) - b_0}{z} &= A(z) - h B(z)
\end{align}
Using the initial values this results in:
\begin{align}
A(z) &= \frac{h z}{1 - (h - 1) z - 2 h z^2} \\
B(z) &= \frac{1 - z}{1 - (h - 1) z - 2 h z^2}
\end{align}
What is left is to express $A(z)$ and $B(z)$ as partial fractions (this will be quite ugly), and read off the coefficients from the resulting geometric series.
