I am looking to finding the solution to indefinite integrals \begin{align} \int z J_m(a z) Y_m(a z)\, \text{d}z, \end{align} and \begin{align} \int z I_m(a z) K_m(a z)\, \text{d}z \end{align}
I know that there are solutions to these when the arguments of these Bessel functions are different, such as the solution posted in Integral of Product of Bessel functions of different kind. Those solutions, though, are singular when the arguments are the same, as is the case I am asking for here.
Are there any ways to get a solution for these integrals? We can assume that $m$ is an integer.
We begin with the indefinite integral
$$\begin{align} \int z Z_m(\alpha z)B_m(\beta z)\,dz=\frac{ \beta z Z_{m}(\alpha z)B_{m-1}(\beta z)-\alpha z Z_ {m-1}(\alpha z)B_m(\beta z)}{\alpha^{2}-\beta^{2}}\tag1 \end{align} $$
for $\alpha \ne \beta$ where $Z_m$ and $B_m$ are arbitrary Bessel functions.
While the left-hand side $(1)$ is well-behaved for $\alpha=\beta$, the right-hand side is not of indeterminate form. The apparent discrepancy is reconciled by remembering that $(1)$ is an indefinite integral. Addition of an appropriate constant with respect to $z$ to the right-hand side of $(1)$ results in an indeterminate form.
Let $V$ be defined by the limit
$$\begin{align} V&=\lim_{\alpha\to\beta}\left(\beta z Z_{m}(\alpha z)B_{m-1}(\beta z)-\alpha z Z_ {m-1}(\alpha z)B_m(\beta z)\right)\\\\ &=\beta z \left(Z_{m}(\beta z)B_{m-1}(\beta z)- Z_{m-1}(\beta z)B_m(\beta z)\right)\tag2 \end{align}$$
The following Wronskians hold:
$$\begin{align} J_{m}(\beta z)Y_{m-1}(\beta z)- J_{m-1}(\beta z)Y_m(\beta z)&=2/(\pi \beta z)\tag{3a}\\\\ H_{m}^{(1)}(\beta z)H_{m-1}^{(2)}(\beta z)- H_{m-1}^{(1)}(\beta z)H_m^{(2)}(\beta z)&=-4i/(\pi \beta z)\tag{3b} \end{align}$$
Given $(3a)$ and $(3b)$, it is evident that $V$, as given on the right-hand side of $(2)$, is independent of $z$. Hence, we subtract $V$ from the numerator of $(1)$ and obtain
$$\begin{align} \int z Z_m(\alpha z)B_m(\beta z)\,dz=\frac{ \beta z Z_{m}(\alpha z)B_{m-1}(\beta z)-\alpha z Z_ {m-1}(\alpha z)B_m(\beta z)-V}{\alpha^{2}-\beta^{2}}\tag4 \end{align} $$
which is clearly of indeterminate form.
It is left as an exercise for the reader to evaluate the limit as $\alpha\to \beta$ of $(4)$ using L'Hospital's Rule to obtain the result
$$\begin{align} \int z Z_m(\alpha z)B_m(\beta z)\,dz&=\frac12 z^2\left(Z_m(\beta z)B_m(\beta z)-Z_{m-1}(\beta z)B_{m-1}(\beta z)\right)\\\\&-\frac{mz}{2\beta}\left(Z_m(\beta z)B_{m-1}(\beta z)+Z_{m-1}(\beta z)B_{m}(\beta z)\right) \end{align} $$