I have an integral function of the form $$ F_n(x) = \int_0^x du \left( K_n(u)L_n(x) - K_n(x)L_n(u) \right)u^p K_n(u) $$ where $$ K_n(x) = x^{(m+1)/2}e^{-x/2}M(a_n,m+1,x), \\ L_n(x) = x^{(m+1)/2}e^{-x/2}U(a_n,m+1,x), $$ $M(a_n,m+1,x)$ and $U(a_n,m+1,x)$ are the confluent hypergeometric functions, $a_n$ is real, $m$ and $p$ are non-negative integers. Hence $F_n(x)$ is real.
I need to evaluate $F_n(x)$ over a large range starting from $x=0$, where $F_n(0)=0$, until a point well in the asymptotic region in which the integrand grows exponentially due to the asymptotic behavior of $K_n(x)$. So, the numerical quadrature I am looking for must work with arbitrary integral limits. In the past, I had used Gauss-Lobatto quadrature in combination with a finite element method but that was for integrating an integrand that decays exponentially.