There are two problems here that need to be solved. The first is the delay of the measured step response, i.e. how many samples from the beginning you need to discard. This is critical in your example because you chose a second-order system to model the data, and such a system is not very flexible in adding or subtracting delay. Delay could be added by adding more zeros to the system, i.e. by increasing the order of the numerator. If you don't want that then you will need to optimize the delay of the step response such that it can be optimally approximated by a second-order system.
The other problem is the fact that the designed filter doesn't reach the final value of the measured step response. The reason is that the impulse response required by Prony's method is derived by computing the difference of adjacent step response samples without taking into account that the step response was truncated. Since several samples of the measured step response are discarded, you also discard the same number of samples of the impulse response, and if you were to compute a step response from that truncated impulse response (by computing the cumulative sum), you would get a step response that is different from the measured one (and that would be close to the one implemented by the designed filter). The solution to this problem is very straightforward. Just truncate the step response and add a value of zero at the beginning. Then compute the corresponding impulse response by taking differences, and the first value of the impulse response will equal the first (non-zero) value of the step response, i.e., it will represent the accumulated past of the original impulse response (before truncation).
I did that and found as optimal starting value for the step response an index of $530$ for a filter order of $2$. The figure below shows that the step response of the designed filter closely approximates the measured step response:
