how to determine rationalfit npoles

I'm trying to figure out how to specify the npoles parameter in rationalfit call... For some fitting the result varies a lot with different npoles value, even for those whose error are below -40dB. So which is the correct value to use? How can I determine that for any given s parameter? Thanks.

To make my point, I copy and make some minor modification to the code on https://www.mathworks.com/help/rf/examples/modeling-a-high-speed-backplane-part-3-4-port-s-parameters-to-differential-tdr-and-tdt.html?searchHighlight=tdr&s_tid=doc_srchtitle

filename = 'default.s4p';
backplane = sparameters(filename);
data = backplane.Parameters;
freq = backplane.Frequencies;
z0 = backplane.Impedance;

wstate = warning('off','rf:rationalfit:ErrorToleranceNotMet');
s11 = rfparam(backplane,1,1);
Vin = 1;
tdrfreqdata = Vin*(s11+1)/2;
[tdrfit p] = rationalfit(freq,tdrfreqdata,'NPoles',1500);
[tdrfit2 p2]  = rationalfit(freq,tdrfreqdata,'NPoles',500);
[tdrfit3 p3]  = rationalfit(freq,tdrfreqdata);
[tdrfit4 p4]  = rationalfit(freq,tdrfreqdata,'NPoles',[0 1500]);
warning(wstate)
Ts = 5e-12;
N = 5000; % number of samples
Trise = 5e-11; % Define a step signal
[Vtdr,tdrT] = stepresp(tdrfit,Ts,N,Trise);
[Vtdr2,tdrT2] = stepresp(tdrfit2,Ts,N,Trise);
[Vtdr3,tdrT3] = stepresp(tdrfit3,Ts,N,Trise);
[Vtdr4,tdrT4] = stepresp(tdrfit4,Ts,N,Trise);

figure
plot(tdrT*1e9,Vtdr,'r','LineWidth',2)
hold on
plot(tdrT*1e9,Vtdr2,'b','LineWidth',2)
plot(tdrT*1e9,Vtdr3,'g','LineWidth',2)
plot(tdrT*1e9,Vtdr4,'c','LineWidth',2)
ylabel('Differential TDR (V)')
xlabel('Time (ns)')
legend('Calculated TDR')

Here is the plot result:

and the fit and error are: p, p2, p3, p4 -54.752132965180394, -53.887304905034920, -26.103337300095674, -41.321671376712070

tdrfit =

   rfmodel.rational with properties:

        A: [1496x1 double]
        C: [1496x1 double]
        D: 0
    Delay: 0
     Name: 'Rational Function'

tdrfit2 =

   rfmodel.rational with properties:

        A: [500x1 double]
        C: [500x1 double]
        D: 0
    Delay: 0
     Name: 'Rational Function'

tdrfit3 =

   rfmodel.rational with properties:

        A: [48x1 double]
        C: [48x1 double]
        D: 0
    Delay: 0
     Name: 'Rational Function'

tdrfit4 =

   rfmodel.rational with properties:

        A: [125x1 double]
        C: [125x1 double]
        D: 0
    Delay: 0
     Name: 'Rational Function'