"plot" command for continous time signal in matlab?uses interpolation? We know digital computer cannot deal with continuous time signals directly. But then how MATLAB "plot" command works?  

Does that uses interpolation technique to show an apparently continous time graph by connecting the various discrete data points?

Till now, i am able to understand is that MATLAB only handles discrete time data, MATLAB models continous time signals using discrete time signals but with very very small sampling interval,so it appears to be continuous while in reality it is discrete and then connects those discrete points using interpolation 

I have also attached a snapshot of wikipedia where a graph of polynomial interpolation is shown 

I have also attached a snapshot of chapter4 first page of book signal processing first, where interpolation discussion is mentioned
By default MATLAB uses linear interpolation when creating line plots, which means it simply draws a line from each point to each point, unless there are more points than pixels in which case each point or group of points within each pixel would represent each pixel shown. 

Compare the following to see this:

This plots a line from x,y=1,2 to x,y= 5,7

plot([1 5], [2 7])


The above is the same as the following where a line is explicitly declared:

plot([1 5], [2 7], '-')


While in comparison the following will plot a “o” at each sample:

plot([1 5], [2 7], 'o') 


Update: To interpolate more smoothly between samples such as polynomial interpolation, you would need to create the interpolated samples and then plot the higher sampled result. The simplest way to do this in MATLAB is to use the resample command where resample(f,M,N) will interpolate the waveform $f$ by $M/N$ samples. Another approach demonstrated below is interpolating using the FFT. Although DSP.SE is not a MATLAB coding site, the following does demonstrate well how to properly do a time domain interpolation using a zero-padded IFFT which is of general signal processing interest. Zero padding in one domain results in interpolation in the other domain, however the technique is to "zero-pad" from the center of the frequency domain waveform with index k from 0 to N-1, such that the high frequency component of the FFT is extended. 

N=6;
M=1000;
x= [0:N-1];
y = sin(2*pi*x/N);
fy = fft(y);
padfy = [fy(1:floor(N/2)) zeros(1, M-N) fy(floor(N/2+1:end))];
figure;
plot(N*[1:M]/M, M/N*real(ifft(padfy)));