Determining Frequency Difference/Drift or Phase Difference/Drift

Hello, The problem that I am trying to tackle is determining frequency/phase difference between two oscillators over time. My two oscillators are GPS locked oscillators (Trimble Mini-T), which each output a 10MHz signal. Looking at both signals on an oscilloscope (Tek MSO2024) I can see that there is indeed some drift between the two over time. They will move in and (slightly) out of phase over periods of time, but the changes are slow.

I have logged data (only 10,000 samples, if you can help me figure out how to take more, that'd be great!) using tmtool in the instrument control toolbox. I took essentially three data sets: some arbitrary phase difference between them, completely inphase, almost completely out of phase.

I have taken these data sets, calculated the FFT of each one, and plotted the results. In each case I get FFTs which indicate that the frequency is precisely 10MHz, for both oscillators, no matter what I do.

Any suggestions of how I could make this measurement? Or perhaps I'm doing something wrong? I will post my code below. I'd like to take longer datasets, but I cannot figure out how to change the "MaxNumberPoint" variable to 100,000. I've tried changing everything in the driver file and I can get my time samples (length(X)) up to 100,000, but never my actual data (X).

At any rate, any help would be much appreciated! Thanks!

Code: if true % code

close all; clear all; clc;
%%Load Captured data from oscilloscope
a = pwd;
d1 = load([a '\data\dataSet_1.mat']);
d2 = load([a '\data\dataSet_2_HandOnOscillator.mat']);
d3 = load([a '\data\dataSet_3_Overcorrection_InPhase.mat']);

Y1 = {d1.Y1, d2.Y1, d3.Y1}; %Channel 1 (Oscillator A)
Y2 = {d1.Y2, d2.Y2, d3.Y2}; %Channel 2 (Oscillator B)

t = d1.t1; %Only need one time set, they are all identical

%Define sampling frequency
Fs = 1e9; %Tek MSO2024 samples at 1GS/s

%%Generate FFT of sinusoids
%generate zero-pad data, making length a power of two

for n = 1:length(Y1)
    nfft = 2^(nextpow2(length(Y1{n})));

      %Take FFT
      fft_y1{n} = fft(Y1{n}, nfft);
      fft_y2{n} = fft(Y2{n}, nfft);

      %Take the positive portion of the FFT
      [fft_y1_pos{n},freqRange{n}] = positiveFFT(Y1{n},Fs);
      [fft_y2_pos{n},freqRange{n}] = positiveFFT(Y2{n},Fs);

      %Plot the FFT
      subplot(3,1,n)
      plot(freqRange{n}/(1e6), abs(fft_y1_pos{n}));
      hold on;
      plot(freqRange{n}/(1e6), abs(fft_y2_pos{n}), 'r');
      grid on;
      xlabel('Freq (MHz)');
      xlim([0 20]);
      hold off;
  end
 end