I am trying to follow a research paper in physics where the authors end up with the following differential equation: $$ D^2(f)=0, \qquad D=\left(\frac{\rm d^2}{{\rm d}t^2} + \frac1r \frac{\rm d}{{\rm d}r} - \frac{1}{r^2}-k^2\right) $$ where the operator $D$ is applied twice, and with $k$ an unknown real number. The article states that using the "Wronskian method" and discarding solutions singular at $r=0$, one can find the following general solution: $$ f(r)=A I_1(kr) + B \left(K_1(kr)\int I_1^2(kr')dr'-I_1(kr)\int I_1(kr')K_1(kr)r'dr'\right) $$ with $A$ and $B$ constant,$I_1$ and $K_1$ the first-order Bessel function of the first and second kind. No more details are given and I'm trying to understand how to do this. I also don't really care about discarding singular solutions at $r=0$.
I can easily find two solutions of this differential equation. Functions satisfying $D(f)=0$, a modified Bessel equation, are also solutions of the initial equation, so we have the two independent family of solutions: $f(r)=\alpha I_1(kr) + \beta K_1(kr)$. Then I assume some tricks based on the Wronskian can be used to find the other two fundamental solutions ...
I looked up examples on how to use the Wronskian but I could't understand how this could be applied here for a 4th order ODE where only 2 solutions are known. From Abel's idenity the Wronskian is $W(r)=k/r^2$, $k$ a constant. I could compute the Wronskian with the determinant definition, but this would give a very ugly expression with 2 unknown functions and their derivatives, which doesn't seem to help ? The resources and exemples I could find always consider an $n$-th order equation for which we already know $n-1$ fundamental solutions, so I'm not sure what to do when two are missing like here.
With $L$ a second order differential operator, the task is to solve $L[L[y]]=0$. Given any solution $y$, then $z=L[y]$ is a solution to $L[z]=0$.
Let $y_{1,2}$ be basis solutions to $L$. Then $z=αy_1+βy_2$. Now it remains to solve $L[y]=z$. Use variation of parameters $$ \begin{bmatrix}y_1&y_2\\y_1'&y_2'\end{bmatrix} \begin{bmatrix}c_1'\\c_2'\end{bmatrix} = \begin{bmatrix}0\\αy_1+βy_2\end{bmatrix} \\ \implies W[y_1,y_2]\begin{bmatrix}c_1'\\c_2'\end{bmatrix} =\begin{bmatrix}y_2'&-y_2\\-y_1'&y_1\end{bmatrix} \begin{bmatrix}0\\αy_1+βy_2\end{bmatrix} $$
This gives a solution formula of $$ y=d_1y_1+d_2y_2-y_1\int\frac{y_2(αy_1+βy_2)}{W[y_1,y_2]}dr+y_2\int\frac{y_1(αy_1+βy_2)}{W[y_1,y_2]}dr $$
In the situation with the singularity the solution $y_2$ that diverges at $r=0$ is discarded, $β=d_2=0$. The Wronskian satisfies $$ W[y_1,y_2](r)\sim \exp\left(-\int_1^r\frac1{r'}\,dr'\right)=\frac1r. $$ This means that there should be a factor $r'$ in both integrands.