The performance of the method is not overwhelming. The bond lengths are overestimated by quite a bit. As cited by Christe, observed is $r_\mathrm{obs}(\ce{I-F_{ax}})=178.1~\mathrm{pm}$ and $r_\mathrm{obs}(\ce{I-F_{eq}})=185.7~\mathrm{pm}$. I calculated $r_\mathrm{BP86}(\ce{I-F_{ax}})=190.4~\mathrm{pm}$ and $r_\mathrm{BP86}(\ce{I-F_{eq}})=195.6~\mathrm{pm}$. Other methods, that I could still reasonably perform did not do much better. The MP2 result is $r_\mathrm{MP2}(\ce{I-F_{ax}})=186.6~\mathrm{pm}$ and $r_\mathrm{MP2}(\ce{I-F_{eq}})=191.9~\mathrm{pm}$.
For the purpose of the analysis the DFT method should be enough. (I had a look at bthethe other results, too; they don't differ significantly.)
I would usually start with looking at the orbitals, but in this case it is more appropriate to look at the electron density itself. Better yet, we look at the Laplacian of the electron density. The depicted QTAIM plots are the Laplacian distribution within the molecule. Solid blue lines indicate charge depletion $\nabla^2\rho>0$, dashed red lines indicate charge accumulation $\nabla^2\rho<0$, red spheres are bond critical points, solid black lines are bond paths, dark red lines are zero flux surfaces (these separate the atoms). See an earlier answer of me for some more information on QTAIM.
We have two bond critical points to consider, one for the axial and one for the equatorial bond. We find low values of electron density at both; $\rho(\ce{I-F_{eq}})=0.125~\mathrm{a.\!u.}$; $\rho(\ce{I-F_{ax}})=0.137~\mathrm{a.\!u.}$. The Laplacian at both points is positive; $\nabla^2\rho(\ce{I-F_{eq}})=0.33~\mathrm{a.\!u.}$; $\nabla^2\rho(\ce{I-F_{ax}})=0.42~\mathrm{a.\!u.}$.
These values indicate predominantly ionic bonding. The graphs nicely show that the atoms appear as separate entities.
With that background it seems wrong arguing the bonding situation by forming covalent bonds, like they are assumed in MO or VB theory.
However, as wrong as it is to neglect ionic bonding, as wrong is it is equally wrong to neglect covalent bonds.
The ionic character of the bonds can also be seen in the charge distribution. While the iodine has a high positive charge $q(\ce{I})=3.24$, the fluorines have a negative charge $q(\ce{F})=-0.46$. These are AIM charges calculated with the MultiWFN program package.
The Analysis with natural bond orbitals give a similar result, $q(\ce{I})=3.20$, $q(\ce{F})=-0.46$.
Looking at the Wiberg bond indices matrix formed from natural atomic orbitals, we find a bond index of 0.54$0.54$ for the equatorial bonds and a slightly higher value of 0.56$0.56$ for the axial bonds.