I'd like to represent the following recursive integral equation to evaluate/graph (for $n\leq3$):
$$K_{i,n}\left(x\right)=\int_{x}^{\infty}K_{i,n-1}\left(y\right)dy$$
where
$$K_{i,0}\left(x\right)=K_{0}\left(x\right)=\int_{0}^{\infty}\frac{\exp\left(-\sqrt{\tau^{2}+x^{2}}\right)}{\sqrt{\tau^{2}+x^{2}}}d\tau$$
The evaluation can be analytic or numeric to serve my purposes. I've tried this:
Kin[x_, 0] := NIntegrate[Exp[-Sqrt[τ^2 + x^2]]/Sqrt[τ^2 + x^2], {τ, 0, ∞}]
Kin[x_, n_] := NIntegrate[Kin[y, n-1], {y, x, ∞}]
which gives:
In[48]:= Kin[1.0, 0]
Out[48]= 0.421024
but also
In[49]:= Kin[1.0, 1]
During evaluation of In[49]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>
During evaluation of In[49]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in y near {y} = {8.16907*10^224}. NIntegrate obtained 3.613500815477093`15.954589770191005*^27949 and 3.613500815477093`15.954589770191005*^27949 for the integral and error estimates. >>
Out[49]= 3.613500815477093*10^27949
These are Bickley-Nayler functions and have definitely been evaluated numerically. The plots from those evaluations look sane and don't indicate that divergence should occur. What have I done wrong in my function definition?