How to avoid high frequency phase error in FFT?

I want to plot the phase response of a signal that is composed of sine waves. After taking the FFT of said signal the phase response features a huge amount of phase error as the frequency axis approaches the sample rate (and even before Nyquist). After some probing I've found that this is somewhat mitigated by ensuring the fft is of length that is a power of 2. The amplitude response the becomes less accurate.

I'm not quite sure why this is happening or how to handle it in different scenarios. I would like the phase response to be smooth over the plotted range. Can anyone tell me what is going on here?

Here is a minimum working example, there is also a

clear

fs = 48e3;

%% Length of signal is power of 2
Ns = 2^16;

t = (0:Ns-1)./fs;
f = (0:Ns-1).*(fs/Ns);

y = sin(2*pi*10*t);

figure(1);
clf;

subplot(211)
loglog(f, abs(fft(y)));

hold on;
xlim([1 fs/2]);

subplot(212)
semilogx(f, unwrap(angle(fft(y))));
hold on;
xlim([1 fs/2]);

%% Length of signal is not power of 2
Ns = fs;

t = (0:Ns-1)./fs;
f = (0:Ns-1).*(fs/Ns);

y = sin(2*pi*10*t);

subplot(211)
loglog(f, abs(fft(y)));

subplot(212)
semilogx(f, unwrap(angle(fft(y))));

legend('Ns is power of 2','Ns is not power of 2','Location','Best');