Unfortunately, I expect that this is too late for the response to be useful for you but I did run into the same issue when investigating the 1996 paper by Rotea et al. that these values come from.
It is my understanding that the thd2thc command is significantly outdated, and has been replaced with d2c for discrete to continuous transfer function conversions. That being said, I don't think that there is any issue with the replacement of the old to new version of this MATLAB function and is simply a precision error caused by the number of decimals included in a1-b3 that you posted. Having 5 significant figures sounds sufficient, but because some of these poles and zero locations are so close to the instability break points on the pole-zero maps, the rounding that you used causes the transfer function to just barely fall into the unstable region. You can prove this by taking Rotea's results after the discrete-to-continous conversion and working backwards to obtain discrete values for comparison...
% NOTATION: ---------------------------------------------------------
% H_ = Transfer Functions originating with the b/a coefficients
% h_ = Computed TF starting with the published Continous results
% _d = Discrete-Time TF
% _c = Continous-Time TF
% _zpk = Zero-pole-Gain TF Format (similar to form in Eq 16 and 17 of
% Rotea et al. paper)
% -------------------------------------------------------------------
% Set the sampling frequency and period for the discrete TF
fs = 50e6; % 50 MHz Sampling Frequency
Ts = 1/fs; % associated sampling period
% Input the coefficients found for equation (11) in Rotea et al.
b1 = -9.2542e-3; b2 = 1.9492e-2; b3 = -1.0238e-2; % Numerators
a1 = -2.9964e0; a2 = 2.9933e0; a3 = -9.9694e-1; % Denominators
% Format num/den form with k = 0 terms
num = [0 b1 b2 b3];
den = [1 a1 a2 a3];
% Generate the discrete TF in digital signal processing (DSP) form
% This would be the same as "tf(num,den,'Variable','z^-1')"
Hd = filt(num,den,Ts);
% Convert this discrete-time transfer function to Continuous-time TF
Hc = d2c(Hd);
% Also convert to zero-pole-gain format
[z,p,k] = tf2zp(Hc.num{1},Hc.den{1});
Hzpk = zpk(z,p,k);
% Plot the step response and bode plots
figure('Name','Step Response from Rotea b/a coefficients')
step(Hc,'r',Hd,'b--',Hzpk,'g-.')
figure('Name','Bode Plots from Rotea b/a coefficients')
bode(Hc,'r',Hd,'b--',Hzpk,'g-.')
% -------------------------------------------------------------------------
% The above results appear to be unstable so work back from the published
% continous-time transfer function presented in Eq 16 of Rotea. It is a
% continuous-time TF presented in zero-pole-gain format
K = -4.8810e5; % Gain
Z = [-1.7152e4; 5.0618e6]; % Zeros of Continuous TF
P = [-3.2607e4; -6.0270e4-1.17e6j; -6.0270e4+1.17e6j]; % Poles of TF
% Generate zero-pole-gain TF model and work backwards
hzpk = zpk(Z,P,K);
% Convert to continuous-time TF
[bnum,aden] = zp2tf(Z,P,K);
hc = tf(bnum,aden);
% Finally Convert to Discrete-time TF
hd = c2d(hc,Ts);
% Plot these with a step response and bode plots as well
figure('Name','Step Response after working backwards')
step(hc,'r',hd,'b--',hzpk,'g-.')
figure('Name','Bode Plots after working backwards')
bode(hc,'r',hd,'b--',hzpk,'g-.')
% -------------------------------------------------------------------------
% Plot the pole-zero maps of the continous and discrete TFs
figure('Name','Pole-Zero Map of Continuous-Time TFs')
pzmap(Hc,'r',hc,'b')
legend('TF using d2c','Published TF after d2c')
figure('Name','Pole-Zero Map of Discrete-Time TFs')
pzmap(Hd,'r',hd,'b')
xlim([0.98 1.03]) % Zoom in to show how close to stable original is
legend('Original Discrete-Time TF','TF working backwards using c2d')
% -------------------------------------------------------------------------
% Print out the transfer functions
brkline = '==============================================================';
abs = '***** TF from plugging in b and a coefficients into eq (11)';
rev = '***** TF working backwards from zero-pole-gain TF in eq (16)';
% Discrete TFs
fprintf('%s\nDISCRETE-TIME TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
fprintf('%s%s%s%s',abs,evalc('Hd'),rev,evalc('hd'))
% Continuous-time TFs
fprintf('\n%s\nCONTINUOUS-TIME TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
fprintf('%s%s%s%s',abs,evalc('Hc'),rev,evalc('hc'))
% Zero-Pole-Gain TFs
fprintf('\n%s\nZERO-POLE-GAIN TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
fprintf('%s%s%s%s\n',abs,evalc('Hzpk'),rev,evalc('hzpk'))
% =======================================================================
% What seems peculiar at this point is that the printed versions for the
% continuous-time TF and the zpk-model TF differ whether we worked forwards
% or backwards, but the discrete-time TF appears to have the same
% coefficients whether we moved forwards or backwards. But investigating a
% little further, it can be seen that this is just due to the rounding on
% the values and this small rounding is what makes all of the difference
% Print the values of in the default rounded format
% First the numerators (top: Forward method. Bottom: Backwards method)
disp([num; hd.num{1}])
% Then the Denominators
disp([den; hd.den{1}])
% The two methods appear to match, but lets print more decimals
format long
% Again the numerator first
disp([num; hd.num{1}])
% Then the Denominators
disp([den; hd.den{1}])
So to sum it up... there is nothing wrong with d2c in comparison to thd2thc but for this specific example, small artifacts from rounding precision makes the original Discrete-time TF unstable and therefore the discrete-to-continuous conversion gives results that don't match what is expected.
