why is the plot of the PSD of this autocorrelation function oscillating?
You are right that the DFT of a real and even sequence is real and even. The problem is that your sequence is not even, at least not the way you programmed it. The fft command doesn't know how you defined the index m. By definition, the first data element in the vector x is associated with the index $n=0$:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$

(where $N$ is the DFT length).

Your definition of the sequence x results in a complex spectrum, and the plot command plots the real part of it.

There are two solutions, one of which you already found: either you take the absolute value of the spectrum, or you define the data in such a way that they are indeed even. Due to the implicit periodicity of the DFT, this means that you append the left half of your data (for negative time indicies) at the end of your data vector. Note that the zero padding must now be done in the middle of your data vector:


x = [x(101:201),zeros(1,Nfft-2*N-1),x(1:100)];
Sx = fft(x,Nfft);
max(abs(imag(Sx))            % 2.4980e-16


The resulting spectrum has a small imaginary part due to numerical inaccuracies, so you always need to use real(fft(x,Nfft)), but it's always good to check if the imaginary part is really as small as it should be (around $10^{-16}$ for double precision floating point format, as used in Matlab).