Starting from physical problem of calculating the vector potential, I come up with following system of differential equations: $$ \begin{cases} \color{red}{R_{zzz}}(1-\lambda^2+z^2)+\color{red}{R_{zz}}(3z)+\color{red}{R_z}(k^2z^2-k^2\lambda^2)+\color{red}{R}(-k^2z)+\color{blue}{E_z}(-k^2z)=\color{orange}{C_z}(k^2z) \\ \color{blue}{E_{zzz}}(1-\lambda^2+z^2)+\color{blue}{E_{zz}}(3z)+\color{blue}{E_z}(1+k^2-k^2\lambda^2)+\color{red}{R_z}(-k^2z)+\color{red}{R}(k^2)=\color{orange}{C_z}(k^2\lambda^2-k^2)\\ \end{cases} $$ where $\color{orange}{C}=\frac{\omega_pq}{ck}\frac{1}{\sqrt{1+z^2}-\lambda z},\quad k\text{ and }z\text{ are free variables}$
Then substitute $\color{red}{R}\to\frac{\omega_pq}{ck}\color{red}{R},\quad \color{blue}{E}\to\frac{\omega_pq}{ck}\color{blue}{E},\quad\color{orange}{C}\to\frac{\omega_pq}{ck}\color{orange}{C}\quad$ and divide both eq. by $\frac{\omega_pq}{ck}$
New $\color{orange}{C}=\frac{1}{\sqrt{1+z^2}-\lambda z}$
Now I look for solution in form of $\color{red}{R}=\sum_{i=0}^{\infty}\color{red}{R^i}\frac{1}{k^i},\quad\color{blue}{E_z}=\sum_{i=0}^{\infty}\color{blue}{E^i}\frac{1}{k^i},\quad$where in case of $\frac{1}{k^i}-$ $i$ means power, in other cases upper index.
It leads to following result:
for odd $i$ holds $\color{red}{R^i}=\color{blue}{E^i}=0$, while for even $i$ there is a reccurent system of differential equations:
$$ \begin{cases} \color{red}{R^i_{zzz}}(1-\lambda^2+z^2)+\color{red}{R^i_{zz}}(3z)+\color{red}{R^{i+2}_z}(z^2-\lambda^2)+\color{red}{R^{i+2}}(-z)+\color{blue}{E^{i+2}_z}(-z)=0\\ \color{blue}{E^i_{zzz}}(1-\lambda^2+z^2)+\color{blue}{E^i_{zz}}(3z)+\color{blue}{E^i_z}(1)+\color{blue}{E^{i+2}_z}(1-\lambda^2)+\color{red}{R^{i+2}_z}(-z)+\color{red}{R^{i+2}}(1)=0\\ \end{cases} $$
With initial conditions $\color{red}{R^0}=0,\quad\color{blue}{E^0_z}=-\color{orange}{C_z}$
It is possible to solve this system in general case:
$$\color{red}{R^{i+2}}=(\frac{1}{\lambda^2})(1+z^2-\lambda^2)^{\frac{1}{2}}\int(1+z^2-\lambda^2)^{-\frac{1}{2}}\left[\color{red}{R^i_{zzz}}(1+z^2-\lambda^2)(-1-z^2+2\lambda^2)+\\+\color{red}{R^i_{zz}}(5z\lambda^2-2z-2z^3)+\color{red}{R^i_{z}}(2z^2+\lambda^2)+\color{red}{R^i}(z)\right]dz+\\+(-\frac{1}{\lambda^2})(1+z^2-\lambda^2)^{\frac{1}{2}}\int(1+z^2-\lambda^2)^{-\frac{1}{2}}\left[\color{red}{R^{i-2}_{zzzzz}}(1+z^2-\lambda^2)^2+\color{red}{R^{i-2}_{zzzz}}(6z)(1+z^2-\lambda^2)+\\+\color{red}{R^{i-2}_{zzz}}(1+10z^2-\lambda^2)+\color{red}{R^{i-2}_{zz}}(3z)\right]dz $$ $\color{blue}{E^i_z}$ is fully defined through $\color{red}{R^i}$: $$ \color{blue}{E^{i+2}_z}=\color{red}{R^{i+2}_z}(z-\frac{\lambda^2}{z})+\color{red}{R^{i+2}}(-1)+\color{red}{R^i_{zzz}}(\frac{1-\lambda^2}{z}+z)+\color{red}{R^i_{zz}}(3) $$
With initial conditions: $\color{red}{R^0}=0,\quad\color{red}{R^2}=(-\frac{1}{\lambda^2})(1+z^2-\lambda^2)^{\frac{1}{2}}\left[\int\color{orange}{C_{zzz}}(z)dz+\int\color{orange}{C_{zz}}(\frac{3z^2}{1+z^2-\lambda^2})dz+\int\color{orange}{C_{z}}(\frac{z}{1+z^2-\lambda^2})dz\right]$
Note: due to physical reasons all constants of integration equal to zero.
And that's where I stuck, having no idea how to get closed form of $\color{red}{R}, \color{blue}{E}$ or at least solve reccurence. Really appreciate if someone can point a way.