DTFT of any finite sequence in matlab I think freqz in a MATLAB toolbox, is the way to obtain DTFT of sequence.
freqz can calculate frequency response of:

H(z)=(Num)/(Den)

We can easily compute the z-transform of any finite sequence x(n) like this:

H(z)=x(0)z^0 + x(1)z^1 + ...

We know in above expression that the Den, is 1.

Recalling that: freeqz(num,den,n) gives the step response in n point.
By x be vector of x(n),

[x1freqz, x1freqzw]=freqz(x,1,3000,'whole');


must gives us  the DTFT.

1)Is it(above statement) correct? what happening if we shift our polynomial?? why?

The second way is to calculate DTFT formula completely, like this:

[X, W]=me_dtft(x1',pi,3000);
figure
title('my')
% plot(W/pi,20*log10(abs(X)));
plot(W/pi,abs(X))
ax = gca;
% ax.YLim = [-40 70];
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

function [X, w]=me_dtft(x,whalfrange, nsample)
    w= linspace(-whalfrange,whalfrange,nsample);
    t=0:1:size(x,2)-1;

    X=zeros(1,size(w,2));
    for i=1:1:size(w,2)
        X(i)=x*exp(-t*1i*w(i))';
    end

end


2)I confused, is the range of parameter t in above code, important?

3)Is this implementation correct? Why?

I think there must be dissonance, since the picture:

Telling us something wrong!
These transform are taken from pure sine wave(its code in in the picture), from right you can see the fft , freqz by the manner explained top and(left) the DTFT as explained earlier.

Edit after "Jason R"'s comment:
Ok, also I removed logarithmic scale, since it make me confuse.
After that, intuitively thy are alike, as you can see in next image, but why they are not exactly the same(refer to last image by logarithmic scale?)?


freqz:

[x1freqz, x1freqzw]=freqz(fliplr(XX'),1,3000,'whole');

figure
title('freqz')
% plot((x1freqzw/pi)-1,20*log10(abs(fftshift(x1freqz))))
plot((x1freqzw/pi)-1,abs(fftshift(x1freqz)))
ax = gca;
% ax.YLim = [-40 70];
ax.XTick = -1:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude')


Sine sample:

Fs=1000; Ts=1/Fs;
time=0:Ts:1;

Freqs=500;
Xs=zeros(length(Freqs),length(time));

for i=1:length(Freqs)
    Xs(i,:)= cos(2*pi*Freqs(i)*time);
end

XX=Xs;
XX=XX./ max(abs(XX));

figure;plot(time, XX); axis(([0 time(end) -1 1]));
xlabel('Time (sec)'); ylabel('Singal Amp.');
title('A sample audio signal');
sound(XX,Fs)
The linspace function as used is going from -fs/2 to +fs/2 for 3000 samples so the +fs/2 is getting counted twice. In contrast the FFT and freqz, which goes from 0 to N-1 do not duplicate the two end points (in this case they go from DC to 1 bin less than fs where fs is the sampling rate). Therefore the sample locations are not exactly the same leading to the difference observable in the two freqz methods.

Further you can get samples of the DTFT by zero padding as another option: fft(x, 3000).

Instead of linspace which will work just fine if you choose the proper start and stop, I like to do:

t = [0: length(x)-1]*1/fs