# A functional recurrence relation with differentiation, closed form

I'm interested in the way of solving the following recurrence relation:

$$a_{n+1}=a_n'+a_1 a_n-b_1 b_n \\ b_{n+1}=b_n'+b_1 a_n+a_1 b_n$$

Where all $$a_n,b_n$$ are functions of $$x$$, and $$'$$ means derivative w.r.t. $$x$$. The initial functions $$a_1$$ and $$b_1$$ are known. Also, assume we have closed forms the derivatives: $$a_1^{(n)}$$ and $$b_1^{(n)}$$ in terms of $$n$$.

I expect there to be a closed form solution for this case, however, I have trouble deriving it.

One idea I had was to use operator notation and matrix notation, then we have:

$$\begin{bmatrix} a_{n+1} \\ b_{n+1} \end{bmatrix}= \begin{bmatrix} D+a_1 & -b_1 \\ D+b_1 & a_1 \end{bmatrix} \begin{bmatrix} a_n \\ b_n \end{bmatrix}$$

Which obviously can be written as:

$$\begin{bmatrix} a_{n+1} \\ b_{n+1} \end{bmatrix}= \begin{bmatrix} D+a_1 & -b_1 \\ D+b_1 & a_1 \end{bmatrix}^n \begin{bmatrix} a_1 \\ b_1 \end{bmatrix}$$

However, the $$n$$th power of a matrix is found through its eigenvalues, and I'm not sure how to find eigenvalues for the operator matrix.

We can also rewrite this in a more clear way:

$$\begin{bmatrix} a_{n+1} \\ b_{n+1} \end{bmatrix}= \left(\begin{bmatrix} 1 & 0 \\ 1 & 0 \end{bmatrix} \frac{d}{dx} +\begin{bmatrix} a_1 & -b_1 \\ b_1 & a_1 \end{bmatrix} \right)^n \begin{bmatrix} a_1 \\ b_1 \end{bmatrix}$$

And use binomial sum? Not sure. How to correctly write binomials sum with a differential operator inside?

Or maybe there's a more simple way to obtain the closed form?

I would prefer to use strictly real methods if possible, thought I kind of doubt it could work without complex numbers.

By letting $$z_n = a_n + i b_n$$ the given equations can be coupled as
$$z_{n+1} = z_n' + z_1 z_n = (D+z_1)z_n = (D+z_1)^n z_1$$
$$z_2=z_1'+z_1^2$$ $$z_3=z_1''+3 z_1 z_1'+z_1^3$$ $$z_4=z_1'''+3 z_1'^2 + 4 z_1 z_1''+6 z_1' z_1^2 + z_1^4$$ and the general form for $$z_n$$ is given by Faà di Bruno's formula.