Here again let me note that exact Nyquist frequency for a pure sine wave should be avoided. Shannon-Nyquist sampling theorem requires that there's no impulse at the Nyquist frequency as a consequence of bandlimitedness, the content at the exact Nyquist frequency is taken to be zero.
Then the following code demonstrates the approximate simulation of an ideal sinc based interpolator applide to (near) critical samples of a pure sine wave. Note that, any finite observation of a signal cannot be bandlimited, so this simulation is not a perfect representation of the true output from an ideal interpolator, nevertheless by chosing signal duration long enough, one can attain an approximately bandlimited signal.
f = 1; % 1 Hz. sine wave...
Fs = 4.2*f; % sampling frequency Fs = 2.2*f ; a bit more than the Nyquist rate.
Td = 25; % duration of observation ultimately determines the spectral resolution.
t = 0:1/Fs:Td; % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end); % get the resulting final duration
L = length(t); % number of samples in the sequence
M = 2^nextpow2(10*L); % DFT / FFT length (for smoother spectral display, not better resolution! )
x = sin(2*pi*f*t); % sinusoidal signal in [0,Td]
%x = x.*hamming(L)'; % hamming window applied for improved spectral display
% Part-II : Approximate a sinc() interpolator :
K = 25; % expansion factor
xe = zeros(1,K*L); % expanded signal
xe(1:K:end) = x;
D = 1024*8;
b = K*fir1(D,1/K); % ideal lowpass filter for interpolation
y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
xlabel('approx simulation of ideal sinc interpolation');
With the result of