I believe this is because the "thd(x)" function calculates the Total harmonic distance in a way different from the one you have manually.
r = thd(x) returns the total harmonic distortion (THD) in dBc of the real-valued sinusoidal signal x. The total harmonic distortion is determined from the fundamental frequency and the first five harmonics using a modified periodogram of the same length as the input signal. The modified periodogram uses a Kaiser window with β = 38.
The following example shows explicitly how to calculate the total harmonic distortion in dBc for a signal consisting of the fundamental and two harmonics. The explicit calculation is checked against the result returned by thd. Notice the additional 10*log() applied when calculating the value manually.
Create a signal sampled at 1 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and two harmonics at 200 and 300 Hz with amplitudes 0.01 and 0.005. Obtain the total harmonic distortion explicitly and using thd.
t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t)+0.01*cos(2*pi*200*t)+0.005*cos(2*pi*300*t);
tharmdist = 10*log10((0.01^2+0.005^2)/2^2)
tharmdist = -45.0515
r = thd(x)
r = -45.0515
References:
Total harmonic distortion - https://www.mathworks.com/help/signal/ref/thd.html#btzx73d_seealso
