interpolating in time domain to improve FFT frequency resolution

Im trying to interpolate time domain data to improve the FFT resolution. Im not sure what Im doing wrong but Im getting a larger amount of bins but not decreasing the frequency spacing between the bins. I think Im doing something wrong or otherwise have made a mistake in theory. It should be simple.


%create a 1MHz sinewave
Fs = 1/1.0000e-09 % Sampling frequency
T = 1/Fs; % Sample time
L = 2500; % Length of signal
t = (0:L-1)*T; % Time vector
x = 15*sin(2*pi*1000000*t)
figure, plot (t, x)


%%FFT before interpolation
Fs = 1/( t(2) - t(1));
L = size(t,2);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
f = Fs/2*linspace(0,1,NFFT/2+1);
test= fft(x ,NFFT)/L;
figure,plot(f,2*abs(test(1:NFFT/2+1)))
ind =find (f>= 1000000,1)
freq_res = (f(2)-f(1))



%interpolate the sinewave by 100
new_t = t(1)t(2)-t(1))/100:max(t);
new= interpn( t, x , new_t, 'bicubic'); %TX

%%FFT before interpolation
Fs2 = 1/( new_t(2) - new_t(1));
L2 = size(new_t,2);
NFFT2 = 2^nextpow2(L2); % Next power of 2 from length of y
f2 = Fs2/2*linspace(0,1,NFFT2/2+1);
test= fft(new ,NFFT2)/L2;
figure,plot(f2,2*abs(test(1:NFFT2/2+1)));
ind2 =find (f2>= 1000000,1)
freq_res = (f2(2)-f2(1))

I was hoping by interpolating the 1MHz sinewave by 1000 I would get a frequency improvement of 1000. Yet, in both, the frequency resolution is approximately the same. Isnt the FFT resolution roughly equal to the sampling frequency divided by the number of samples. I know we are increasing the number of samples, but isn't interpolation similar to sampling at higher frequency?