Solution of an inhomogeneous modified Bessel equation I'm solving the equation

$$x^2y''+xy'-(x^{2}\lambda^{2}+1)y=-C\frac{I_{1}(\lambda x)}{I_{1}(\lambda)}$$

where $I_{\alpha}$ is the modified Bessel function of the first kind.
The complementary function for this problem is
$$y_{cf}=AI_{1}(\lambda x)+BK_{1}(\lambda x),$$
where $K_{\alpha}$ is the modified Bessel function of the second kind.
Using the method of variation of parameters I find that
$$y_{p}=u_{1}I_{1}(\lambda x)+u_{2}K_{1}(\lambda x),$$
where
\begin{align*}
u_{1}&=\frac{C}{I_{1}(\lambda)}\int\frac{I_{1}(\lambda x)K_{1}(\lambda x)}
{W[I_{1}(\lambda x),K_{1}(\lambda x)]}\,\textrm{d}x
=-\frac{C}{I_{1}(\lambda)}\int xI_{1}(\lambda x)K_{1}(\lambda x)\,\textrm{d}x,
\\
u_{2}&=-\frac{C}{I_{1}(\lambda)}\int\frac{I_{1}(\lambda x)I_{1}(\lambda x)}
{W[I_{1}(\lambda x),K_{1}(\lambda x)]}\,\textrm{d}x
=\frac{C}{I_{1}(\lambda)}\int xI_{1}(\lambda x)I_{1}(\lambda x)\,\textrm{d}x.
\end{align*}
Here I have used the fact that
$$W[I_{1}(\lambda x),K_{1}(\lambda x)]=I_{1}(\lambda x)[K_{1}(\lambda x)]'-K_{1}(\lambda x)[I_{1}(\lambda x)]'=-\frac{1}{x}.$$
Can anyone suggest the best way to go about computing these integrals? I've tried 'by parts' (seemed the natural choice) but haven't managed to get anywhere thus far - thanks!
 A: Consider the following differential equation
$$\frac{{{d}^{2}}u}{d{{z}^{2}}}-{{k}^{2}}u-\frac{q\left( q+1 \right)}{{{z}^{2}}}u=0$$ 
Note firstly that if we set $y=\frac{u}{\sqrt{z}}$ then we have
$${{z}^{2}}y''+zy'-\left( {{k}^{2}}{{z}^{2}}+{{\left( q+\tfrac{1}{2} \right)}^{2}} \right)y=0$$
 Then we set $q+1/2=\nu ,\,\,ikz=x$and so obtain
$${{x}^{2}}y''+xy'+\left( {{x}^{2}}-{{\nu }^{2}} \right)y=0$$
 Bessel’s DE.  Hence $y={{Z}_{\nu }}\left( x \right)$ where Z is the cylindrical function of interest and of order $\nu $.  Inverting the transform therefore we see $u={{z}^{1/2}}{{Z}_{q+\tfrac{1}{2}}}\left( ikz \right)$ .  In a similar fashion if we set $u=w{{z}^{-q}}$ we obtain
$$\frac{{{d}^{2}}w}{d{{z}^{2}}}-\frac{2q}{z}\frac{dw}{dz}-{{k}^{2}}w=0$$
Which has a solution $w={{z}^{q+1/2}}{{Z}_{q+\tfrac{1}{2}}}\left( icz \right)$ .  We now get to a result due to Lommel, one which I will simply state.  Replace $w=y/\theta \left( z \right)$ and $z\to \psi \left( z \right)$ and $2q$ by $2\nu -1$ and $k=i$ , then 
$$\begin{align}
  & \frac{{{d}^{2}}y}{d{{z}^{2}}}-\left( \frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta } \right)\frac{dy}{dz} \\ 
 & +\left\{ \left( \frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta } \right)\frac{\theta '}{\theta }-\frac{\theta ''}{\theta }+\psi {{'}^{2}} \right\}y=0 \\ 
\end{align}$$ 
The solution of which is $y=\theta \left( z \right)\psi {{\left( z \right)}^{\nu }}{{Z}_{\nu }}\left( \psi \left( z \right) \right)$ .  Now introduce a new function defined by
$$\frac{\varphi '}{\varphi }=\frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta }$$ 
Using this we eliminate $\theta $ to obtain
$$\begin{align}
  & \frac{{{d}^{2}}y}{d{{z}^{2}}}-\frac{\varphi '}{\varphi }\frac{dy}{dz} \\ 
 & +\left\{ \frac{3}{4}{{\left( \frac{\varphi '}{\varphi } \right)}^{2}}-\frac{1}{2}\frac{\varphi ''}{\varphi }-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\frac{1}{2}\frac{\psi '''}{\psi '}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}y=0 \\ 
\end{align}$$
The solution of which has the form
$$y=\sqrt{\frac{\varphi \psi }{\psi '}}{{Z}_{\nu }}\left( \psi  \right)$$ 
Now let $\varphi =1$ and hence
$$\frac{{{d}^{2}}y}{d{{z}^{2}}}+\left\{ \frac{1}{2}\frac{\psi '''}{\psi '}-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}y=0$$
Now observe that if ${{y}_{\nu }}''+P{{y}_{\nu }}=0$ and ${{y}_{\mu }}''+Q{{y}_{\mu }}=0$ then we have 
$${{y}_{\nu }}{{y}_{\mu }}''-{{y}_{\nu }}''{{y}_{\mu }}=\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}$$
or
$$\frac{d}{dx}\left( {{y}_{\nu }}{{y}_{\mu }}'-{{y}_{\nu }}'{{y}_{\mu }} \right)=\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}$$
or
$$\int{\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}dx}={{y}_{\nu }}{{y}_{\mu }}'-{{y}_{\nu }}'{{y}_{\mu }}$$
Here we have 
$$P,Q=\left\{ \frac{1}{2}\frac{\psi '''}{\psi '}-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}$$ 
We now let $\psi =az,\,\,\nu =\nu $ for P, and $\psi =bz,\,\,\nu =\mu $for Q.  This yields 
$$P=\left( {{a}^{2}}{{z}^{2}}-\nu +\frac{1}{4} \right)\frac{1}{{{z}^{2}}},\,\,Q=\left( {{b}^{2}}{{z}^{2}}-\mu +\frac{1}{4} \right)\frac{1}{{{z}^{2}}}$$
Let${{y}_{\nu }}=\sqrt{z}{{Z}_{\nu }}\left( az \right)$and ${{y}_{\mu }}=\sqrt{z}{{\bar{Z}}_{\mu }}\left( bz \right)$ where the bar on Z here indicates that it is possibly a different cylinder function than that in ${{y}_{\nu }}$.  We have then
$$\begin{align}
  & \int{\left( \left( {{a}^{2}}-{{b}^{2}} \right)z+\left( \mu -\nu  \right)\frac{1}{z} \right){{Z}_{\nu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( bz \right)dz} \\ 
 & =z\left( {{Z}_{\nu }}\left( az \right)\frac{d}{dz}{{{\bar{Z}}}_{\mu }}\left( bz \right)-{{{\bar{Z}}}_{\mu }}\left( bz \right)\frac{d}{dz}{{Z}_{\nu }}\left( az \right) \right) \\ 
\end{align}$$
Now let $\nu =\mu $ and so
$$\int{z{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( bz \right)dz}=\frac{z\left( {{Z}_{\mu }}\left( az \right)\frac{d}{dz}{{{\bar{Z}}}_{\mu }}\left( bz \right)-{{{\bar{Z}}}_{\mu }}\left( bz \right)\frac{d}{dz}{{Z}_{\mu }}\left( az \right) \right)}{\left( {{a}^{2}}-{{b}^{2}} \right)}$$
Using L’Hopital’s rule we have
$$\int{z{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( az \right)dz}=\frac{{{z}^{2}}}{4}\left( 2{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( az \right)-{{Z}_{\mu -1}}\left( az \right){{{\bar{Z}}}_{\mu +1}}\left( az \right)-{{Z}_{\mu +1}}\left( az \right){{{\bar{Z}}}_{\mu -1}}\left( az \right) \right)$$
Therefore, for example
$$\int{z{{I}_{1}}\left( az \right){{I}_{1}}\left( az \right)dz}=\frac{{{z}^{2}}}{2}\left( {{I}_{1}}{{\left( az \right)}^{2}}-{{I}_{0}}\left( az \right){{I}_{2}}\left( az \right) \right)$$
All of this is in Watson's treatise (pg 134-135), although a lot more detail has been added here.
