Let $r,L$ be the radius and length, and let $K$ be the surface current density -- i.e., $K=(N/L)I$ where $N$ is the number of turns. The standard, approximate formula for the energy stored is
$$
U_0 = {1\over 2}\mu_0K^2\pi r^2L.
$$
Let $U_{\rm in},U_{\rm out}$ be the exact expressions for the energy stored inside and outside the solenoid:
$$
U_{\rm in}={1\over 2\mu_0}\int_{\rm interior} B^2 d\tau,
$$
where $d\tau$ is a volume element, and the integral is taken over the inside of the solenoid,
and similarly for $U_{\rm out}$.
I claim that, for $L\gg r$, difference between $U_{\rm in}$ and $U_0$ is approximately
$$
U_{\rm in}-U_0\approx -c_1U_0{r\over L}.
$$
where $c_1$ is some constant of order 1. Here's why. For large values of $L/r$, the magnetic field inside the solenoid will be almost exactly equal to the value you get for an infinite solenoid (i.e., $B=\mu_0K$). The deviations will occur only near the ends, in a region of order $r$. The field in this region near the end cap will not depend on $L$ -- near one end, the solenoid "looks" semi-infinite. So that energy difference should depend on $r$ but not $L$. Moreover, it must be proportional to $\mu_0K^2$. On dimensional grounds, then, it's got to be of the form $\mu_0K^2r^3$ times some dimensionless constant. That's equivalent to the form I have above.
Oh, and the field gets weaker near the ends of the solenoid -- that's why I put in that minus sign.
The same sort of argument applies to the energy outside the solenoid. The energy stored outside the solenoid is
$$
U_{\rm out}\approx c_2U_0{r\over L},
$$
where $c_2$ is some dimensionless constant of order 1.
The fraction of the total energy stored inside the solenoid
$$
{U_{\rm in}\over U_{\rm in}+U_{\rm out}}={1-c_1(r/L)\over 1+(c_2-c_1)(r/L)}\approx 1-c_2(r/L)
$$
(the approximation being the leading term in an expansion in $r/L$ as usual). That's the fraction of the actual total energy. The fraction of the nominal total energy is
$$
{U_{\rm in}\over U_0}=1-c_1(r/L).
$$
Either way, the rfaction approaches 1 for large $L$.
Now all that's left is to find $c_1$ and $c_2$. This is messy. The magnetic field of a finite solenoid can be written out exactly in terms of elliptic integrals. Just square them and do the appropriate volume integrals, and you have it.
Actually, the easiest thing to work out from these formulas is the total energy in all space,
$U_{\rm in}+U_{\rm out}$. The reason is that this one can be computed from the vector potential, by integrating ${\bf J}\cdot{\bf A}$, rather than from the field. It turns out that
$$
U_{\rm in}+U_{\rm out}=U_0\frac{\sqrt{\lambda ^2+4} \left(\lambda ^2 K\left(\frac{4}{\lambda
^2+4}\right)-\left(\lambda ^2-4\right) E\left(\frac{4}{\lambda
^2+4}\right)\right)-8}{3 \pi \lambda }
$$
where $\lambda=L/r$. Taylor expanding in $1/\lambda$, we get
$$
U_{\rm in}+U_{\rm out} \approx U_0\left(1-{8\over 3\pi}{r\over L}\right)
$$
for large $L$. In other words,
$$
c_2-c_1=-{8\over 3\pi}=-0.8488.
$$
To get the $c$'s individually, the best thing I could find to do was to numerically integrate $B^2$ in the interior of the cylinder. This results in $c_1\approx 1.50$, which means $c_2\approx 0.65$.