The following recurrent relation is posed
$\begin{cases} b_{i+1} = (1-\frac2N)b_i+2(1-\frac1N)a_i,\\ a_{i+1} = -\frac2Nb_i + (1-\frac2N)a_i \end{cases}$
Now I would like to understand how the author of the paper I referenced below (Section 4.3, page 8 of pdf, page 472 of journal, http://www.complex-systems.com/abstracts/v11_i06_a03.html) solved this problem. I already contacted him, but I received no reply. I understand that he determines the solution by first transforming the recurrent relations to a difference between two steps
$\begin{cases} b_{i+1} - b_i = \frac2Nb_i+2(1-\frac1N)a_i,\\ a_{i+1} - a_i = -\frac2Nb_i + -\frac2Na_i \end{cases},$
then extends it to difference equations
$\begin{cases} \dot{b} \delta = \frac2Nb_i+2(1-\frac1N)a_i,\\ \dot{a} \delta = -\frac2Nb_i + -\frac2Na_i \end{cases}$
for $\delta$ small as required. Now I do not understand what he does in the next step, when he derives
$$\ddot{b} + \frac4{\delta^2 N}b+O(\frac{b}{\delta^2N^2})+\ddot{b}O(\frac1N)+\dot{b}O(\frac1{N\delta})=0$$
Also it seems for me strange to keep the two terms $\frac4{\delta^2 N}b+O(\frac{b}{\delta^2N^2})$, since $\frac4{\delta^2 N}b+O(\frac{b}{\delta^2N^2})\in O(\frac{b}{\delta^2N^2})$
The way I know to solve a linear system of differential equations is to formulate the equations into a linear system, determine eigenvalues and -vectors, and insert them into a linear combination of eigenevectors times exponential functions like it is described in tutorial.math.lamar.edu/Classes/DE/RealEigenvalues.aspx (cannot post link because of reputation system), but I do not have any clue how the author comes to this equation.
Ozhigov, Yuri, Speedup of iterated quantum search by parallel performance, Complex Syst. 11, No.6, 465-486 (1997). ZBL0955.68036.
I am not sure if you need to convert the recurrence system into a differential equation to solve. Define $c_i = \begin{pmatrix}b_i\\a_i\end{pmatrix}$.
Let $\delta = 1-{2\over N}$, and $A = \begin{pmatrix}\delta \ \ \ \ \ \delta+1 \\\delta-1 \ \ \ \ \ \delta \end{pmatrix}$.
Then,
$$c_{i+1}=Ac_i, $$ or $$c_k=A^kc_0.$$
To compute $A^k$, you can try decomposing the matrix using eigenvectors, $v$, such that $(A - \lambda I) v = 0$
It is a bit of algebra, but the eigenvectors and eigenvalues are: $$v_1 = \begin{pmatrix}i\sqrt{\frac{1+\delta}{1-\delta}}\\1\end{pmatrix}, \lambda_1=\delta-i\sqrt{1-\delta^2}$$ $$v_2 = \begin{pmatrix}-i\sqrt{\frac{1+\delta}{1-\delta}}\\1\end{pmatrix}, \lambda_2=\delta+i\sqrt{1-\delta^2}$$
So, $A^k = V\Lambda^k V^{-1}$, where $\Lambda=\begin{pmatrix}\lambda_1\ 0\\0 \ \ \lambda_2\end{pmatrix}$, and $V=[v1, v2]$
Addendum: This simplifies further (again, with more algebra). Let $\phi = \tan ^{-1}\frac{\sqrt{1-\delta^2}}{\delta}$ and $\beta = \sqrt{\frac{1+\delta}{1-\delta}}=\sqrt{N-1}.$
Then, $$c_k=\begin{pmatrix}\cos\phi k \ \ \ \ \ \ \beta\sin \phi k \\ -\beta^{-1}\sin \phi k \ \ \ \ \cos\phi k\end{pmatrix}c_0$$
For $N>>1$, $\beta\sin \phi k \approx 2k\frac{N-1}{N-2}$. As a sanity check, in the limit as $N$ approaches $\infty$, notice the recurrence becomes $a_n = a_0$, and $b_{n+1}=b_n+2a_0$. Thus, $b_n=2na_0+b_0$