The feff document tells us this about the IORDER parameter (links and some formatting added by me):

Order of the approximation used in module GENFMT. feff uses order 2 by default, which is correct to terms of order 1/(pR)^2 and corresponds to 6x6 scattering matrices in the Rehr-Albers formalism. Single scattering is calculated exactly to this order. The 6x6 approximation is accurate to within a few percent in every case we have tried (that is, higher order doesn’t change the result more than a few percent). However M4 shells and higher shells may require increased iorder for coupling the matrix elements. Changing the default values requires some familiarity with the Rehr-Albers paper and the structure of the module GENFMT. To do so, follow the instructions in the feff source code in subroutine setlam. The key iord is passed to setlam for processing. You may need to change the code parameter lamtot if you want to do higher order calculations. For details of the algorithm used by GENFMT, see the paper by J.J. Rehr and R.C. Albers. For the M4 and higher edges, you may receive an error message like: Lambda array overfilled. In that case the calculations should be repeated with IORDER -70202 (10x10 matrices).

To test the effect of changing the iord parameter on EXAFS analysis, I compiled up a copy of the genfmt program with /src/HEADERS/dim.h modified with lamtot=35, mtot=6, and ntot=4. (I am a bit skeptical that I have done this correctly. I can run genfmt to completion with the default values (15,4,2) and get the same results. Either the advice to modify the code for higher order is nonsense, or it is not explained clearly enough for me to follow.)

Caveat: I selected those values based on my understanding of setlam.f. genfmt ran to completion without complaint, so I am hopeful that that was done correctly.

For each material (see the SCF tests document for descriptions of the materials), I computed Feff with self-consistency and the self-consistency radius set to the second shortest value used in the SCF tests. For example, for FeS2, the radius was set to 3.6 Å and, for BaZrO3, the radius was set to 4 Å. I then ran calculations with the iord parameter set to 1, 2, 3, 4, and 10. The code identifies 10 as triggering the “cute” algorithm, which treats collinear paths differently from other multiple scattering paths.

Changes to the iord parameter should only effect multiple scattering paths. Single scattering paths are calculated without that approximation. This is easily tested. Running a sequence of first shell fits with different values if iord does, in fact, result in identical fit results. For example, here are the results for first shell fits to FeS2:

1. A couple of the materials behave pretty much as one might expect. NiO, Bromoadamantane, and BaZrO3 show a significant drop in chi-reduced between iord of 1 and 2, while not showing a statistically significant change in any of the fitting parameters.
2. A few materials – Copper, FeS2, and uranyl – actually show somewhat better chi-reduced for iord of 1. I think this tells us that at iord of 1, the calculation is not converged and that the effect of this on the EXAFS analysis is ill-determined. I think it would be a mistake to claim something like “fitting is better in some cases with iord=1.” Rather, this variability is telling us that iord=1 is a mistake.
3. In most cases, there is very little change in chi-reduced for iord>=2. While there is some variability among the larger iord results for some materials (NiO, for example, varied by a bit more than 1%), it seems that the default of iord=2 is well justified.
4. Perhaps this exercise could be used to approximate the systematic uncertainty contributed by the MS theory to the EXAFS analysis….