Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples
It's your observation interval which creates the main problem.

Your reasoning based on the Nyquist sampling theorem is ok; of course with a pure sine wave at the exact Nyquist frequency you will have troubles and therefore it's wise to relax the sampling frequency (slightly) above that of Nyquist rate, such 2.2 Hz instead of a strict 2 Hz... So this is one problem you will practically have. 

But your main problem that appears on the FFT plot is about the spectral resolution due to short observation interval. Mainly with two samples (or one second observation of a 1 Hz sine wave) your FFT resolution will be limited to even less than one Hz. Please search the site for PSD, spectral resolution, FFT bin frequency to get a better understanding of spectral analysis of practical windowed data.

In order to see sharp frequency peaks (ideally impulses) in the FFT output of your sine wave, you should increase the spectral resolution, which requires you to increase the observation interval.

I have modified and also expanded your code to see the result of ideal sinc based interpolator on the (near) critically sampled data. Note that I've mereley included a digital simulation of ideal sinc based interpolator (not the simulation of analog interpolator) to see that it would indeed reconstruct the pure sinudosidal from its given samples taken at close to Nyquist rate. Note that for the ideal sinc interpolator to work, the original signal should strictly be bandlimited, or at least sufficiently so, which will have a lot of consequences on the success and efficiency of interpolation.  

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);

subplot(3,1,1);
plot(t,x);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');


subplot(3,1,3)
plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');


Below is a plot for the result of interpolation from a set of near critical sampling.



And below is the same simulation with a more relaxed sampling, as you can see the interpolator performs much better for this set of improved samples (consequence of better bandlimitedness)