The formula you're using is wrong -- that's the base case by dimensions alone (i.e. you put together a length and area, and permeability, and you get an inductance), and there is a constant factor due to geometry of the problem.
As it happens, the geometry factor is 1 for a "long solenoid", or a "thin" toroid.
For a finite rod core like this, it takes on other values.
Consider the old Wheeler formula,
$$ L = \frac{r^2 N^2}{9 r + 10 l}$$
for L in µH, turns \$N\$, and dimensions (coil mean radius \$r\$ and on-centers length \$l\$) in inches (a different coefficient applies for SI units!). This formula is quite simple, and isn't too bad of a fit considering; it's within 10% or so for most practical values.
There's also a multilayer formula, see H. A. Wheeler, "Simple Inductance Formulas for Radio Coils," in Proceedings of the Institute of Radio Engineers, vol. 16, no. 10, pp. 1398-1400, Oct. 1928.
Applying relative permeability \$\mu_r\$ scales the value proportionally, but notice this assumes all the air around the coil has been replaced by the permeable material.
But that's not the case for the rod core, for which there's a whole surrounding air space still.
So we need different formulas.
If simple formulas weren't readily apparent back in the day, the alternative was to plot a whole bunch of points (whether hard-won by empirical work, or laborious calculation of complicated formulas) and prepare a graph or nomogram. We have that here:
(I have this from Pressman, Billings, Morey, Switching Power Supply Design, 3rd Ed., Mcgraw-Hill (2009), but their source in turn comes from Micrometals, see: https://www.tme.eu/Document/82fc93ef93d05a020041f9203a13131a/MICROMETALS-Tseries.pdf page 24 (26). We also notice a discrepancy: they give multilayer radius as to the outer layer, but this is erroneous, it should be the average diameter.)
Note that the diagram doesn't extend to very high permeability such as ferrite provides; presumably the family of curves converge to a vertical asymptote, but it sure looks like the large L/D curves could reach a high value \$\mu_\textrm{eff}\$ on their way there.
Note also, the winding covers the whole core; in general, the inductivity (\$A_L\$ in µH/t2 or whatever) varies with the position of a given turn on the core, and turns don't couple perfectly with each other. The average figure, across all turns along the length of the rod, is as given here (or well, after dividing by \$N\$). You will get different results for short windings that don't cover the whole rod.
I don't have a close calculation for your particular materials and geometry, but a deeper lesson is perhaps this: magnetics, and fields problems in general, are often much more complicated to solve than you might expect given the simple materials needed to construct them; as such, either mere crude approximate solutions exist, or you basically have to go and measure the whole thing -- whether numerically in a field simulator, or empirically (the real deal).