Hi,
your calculation of 3.6MOhm is correct, as the differential output resistance is defined as R = dU / dI.
Also your method is ok.
This type of GND referenced current source has the disadvantage that it needs matched resistors, but you might trim one of them (R1) to increase the output resistance to near infinity.
A further disadvantage is the common mode voltage on the OpAmp, which might be included in the calculation.
I did not examine your calculation, but you better should calculate the total differential (for your own exercise

) which should give a more clear mathematical expression for the output impedance over mismatch.
You do not necessarily need a 3458A, (ok it has very versatile DCI ranges, but 7.5 digits only in that mode)
I assume, you did not set up the 3458A correctly, much higher values of output resistance should be possible to measure with that box...
Therefore any 6.5 digit, high impedance DMM as a 34401A will also do the job:
Use it in DCV mode, high impedance, put an appropriate stable shunt resistor in parallel, 2k in this case, yielding max. resolution at 1.000000 V for 500µA, i.e you can determine changes as low as 500pA on the last digit, while switching between 0 Ohm and 10k for a 5V voltage swing (instead of only 0.3V for 600Ohm), in series with the instrument .
Using +/-15V for the OPAMP should give enough volt compliance.
If the current change happens on the last digit only, you would measure 1,000,000 * 10k = 10GOhm of diff. output resistance.
That's geeky enough, isn't it?
Frank