The above will work if the winding is confined to a small fraction of the length, so that leakage around the rod largely avoids the winding itself.
For a distributed winding, significant leakage will span the winding itself, and the effective area will be larger. The field will still be strongest in the middle, but the number of turns acting exactly upon the middle is fewer, and the coupling from more distant turns is less than one. It will have a higher saturation current, but it's hard to say by how much.
In the extreme, for a rod much shorter than the winding, the reduction of magnetic path length is very small indeed, so we might use the solenoid formula to find the flux density inside the winding and assume it saturates when B_air = B_sat. But this wouldn't be very useful as the inductance gain is also tiny (equivalent to the given, that the core length is a small fraction of the total magnetic path length [through air]).
I think I would estimate by starting with the mu_eff estimate,

and assume that mu_eff gives the average reduction in magnetic path length, and so we can assume the same ratio gives the increase in flux density over the air-cored case.
Still very rough, but perhaps good enough for use.
Easy enough to wind a few and measure saturation, at least; make a scale model if necessary.
FWIW, single turns on powdered iron toroids (mu >= 50, say) are about the right ballpark there. Might use them like ferrite beads?
Tim