How to find phase difference between 2 signals, one has an unknown function
I don't fully understand your question, or what your code is supposed to do. However, this line is wrong:

fs =  14*frequency;


You should use the same sampling frequency as was used to acquire the data. Assuming uniform sampling, you can do something like

dt = t(2) - t(1);
fs = 1/dt;


Also, be aware that estimating the actual phase from the FFT is not straigthforward. Here is a quick and dirty example showing one way to do it. Start with this code:

fs = 1000;
t = 0:1/fs:0.1;
f = 100;    % signal frequency
phase = 0;  % signal phase
y = cos(2*pi*f*t+phase);
Y = fftshift(fft(y,2^nextpow2(length(y))));
Ym = abs(Y);
[m, idx] = max(Ym);
Ym(idx) = 0;
[m, idx] = max(Ym);  % index of f
Ya = angle(Y);
ph_e = Ya(idx);
fprintf('Actual phase is: %1.5f\n', phase);
fprintf('Estimated phase is: %1.5f\n', ph_e);


If you run this code you'll get:

Actual phase is: 0.00000
Estimated phase is: -0.49087


The phase estimate depends on many factors, including the start and end times of your time vector. The thing is that, in this particular setup, the FFT interprets a phase of zero as -0.49087. Then, modify the third to last line in the script:

ph_e = Ya(idx) + 0.49087;


Now the script will properly identify the phase. Try it by setting the phase in line 4 to, say, pi/3, then run the script:

Actual phase is: 1.04720
Estimated phase is: 1.03347


Which is pretty good accuracy.

I hope this helps you to get going.