The original question talked about a discrepancy between the result obtained by calculating the flux directly and using the definition $L= \Phi/I$. The confusion arises because of a concept known as "flux linkage". When you calculate the flux enclosed by the region of unit length between r and r+dr, you calculated an expression for flux which you integrated to get the total flux. However, the entire flux calculated by you is not "linked" to this area since the current enclosed by the contour of this radius is a fraction $\pi r^2/R^2$ of the total current. Thus the linked flux is $d\Phi= \dfrac{\mu_0 I r dr}{2\pi R^2} \dfrac{\pi r^2}{\pi R^2}$. If you integrate this expression, you would get the correct result. $$\Phi = \dfrac{\mu_0 I}{2\pi R^4}\int_0^\infty r^3 dr= \mu_0 I/8\pi$$$$\Phi = \dfrac{\mu_0 I}{2\pi R^4}\int_0^R r^3 dr= \mu_0 I/8\pi$$ "Flux linkage' is not a very easy concept but consider what happens when you have N turns of the wire through which the same flux passes. In order to calculate emf using Faraday's law, you will need N times the flux to get the correct emf.
Changed upper limit of integral to R, from inf, which it was. Int (r^3) to infinity is undefined.
