The coherence equals one in your example because you estimated the spectra Pxx, Pyy, and Pxy using a single segment. If you use multiple segments, the coherence esitmate will have be between 0 and 1. I will explain this more below, but first:
The time vector in your data is sampled at a non-uniform rate. See plot of delta t versus sample number below.
The plot shows that there are handful of points sampled fast, then about 60% of the points are sampled relatively slowly, then the last 40% (approximately) of points are sampled about 17 times faster (dt is 17 times smaller) than the first 60%. Therefore you should resample to a uniform rate, with interp1(), before estimating spectra. Or maybe the time vector in your data file is wrong.
Let's pretend the data is sampled at a uniform rate, with dt=1, i.e. fs=1/dt=1.
Estimate coherence the way you did, with a single segment that includes all the data.
Coh1=abs(Pxy1).^2./(Pxx1.*Pyy1);
Next: estimate coherence with segments that are one-eighth of full length (or just a bit shorter than 1/8, so that 8 segments will fit), plus 7 segments of equal length that overlap the first 8 by half on each side. Thus there will be 15 total segments, each with length <= 1/8 of the full length sequence. This approach uses all or almost all the data, and reduces the variance in the power spectral estimate, compared to the estimate obtained with a single full-length segment.
if mod(Nseg,2)==1,Nseg=Nseg-1; end
Pxx8=cpsd(x,x,Nseg,Nseg/2,Nseg);
Pyy8=cpsd(y,y,Nseg,Nseg/2,Nseg);
Pxy8=cpsd(x,y,Nseg,Nseg/2,Nseg);
Coh8=abs(Pxy8).^2./(Pxx8.*Pyy8);
Plot results
plot(f1,Coh1,'-r.',f8,Coh8,'-b.')
grid on; xlabel('Frequency'); ylabel('\gamma^2'); title('Coherence Estimates')
legend('1 segment','15 segments');
The code above produces the figure below:
The figure shows that when one estimates coherence with one segment, the coherence equals 1 at all frequencies. But if one uses multiple segments, the estimate of coherence is more realistic, and is between 0 and 1. Frequency goes from 0 to 0.5 in the plot above, because we assumed sampling rate fs=1. Thefore the coherence spectra go up to the Nyquist frequency.
Computing coherence with a single segment is kind of like like estimating Pearson correlation (r-squared) for the best-fit straight line, for a data set with just two points. The line fits the data perfectly when there are only two points, so r^2=1.000. Not a perfect analogy, but close.