Finding phase difference using cross correlation

Hi, I am finding phase difference (time lag) between two sine waves. My sampling frequency is 500K and I am interested in frequencies 5k, 6k, 7k.......50k. I am using this simple code:

 fs=500000; 
 t = [0:1/fs:0.001];
 A = 1;
 f =5000;
 y = A*sin(2*pi*f*t);
 figure(1)
 plot(t,y)
 y1 = A*sin(2*pi*f*(t-0.00005));
 figure(2)
 plot(t,y1)

 x = xcorr(y,y1,'coeff');
 tx = [-(length(y)-1):length(y1)-1]*(1/fs);
 [mx,ix] = max(x);
 lag = tx(ix);

I am getting the accurate value of lag for lower frequencies like 5k to 20k. But as frequency increases, the value of lag is not accurate. When I increased my sampling frequency to 5000K (from 500k), I am getting the accurate value of lag for frequencies 5k to 40k. Finally, when my sampling frequency is 7000k, I am getting accurate lag for all the frequencies (5k to 50k).

Is there any way to calculate accurate lag at the low sampling frequency (around 500k)? Is correlation technique sufficient for this purpose? If not, then what can be the other techniques?

[Ernst Theussl]

In your code you've defined signal with the phase shift via

y1 = A*sin(2*pi*f*(t-0.00005));

But I think it should be:

y1 = A*sin(2*pi*f*t-0.00005);