Adapt fft to goertzel algorithm

What Im trying to do is using FFT and Goertzel algorithms to find the phase difference between two sin signals. I believe that my code run just fine with FFT. It can output the correct phase difference value, but the Goertzel algorithm cannot output the expected result. From what I understand, the difference between Goertzel algorithm and FFTis that the Goertzel algorithm has a region of interest. I don't know why when I am doing the Goertzel algorithm, it only has one value for each sin function. The result screenshot and code are shown below.

clear;

clc;

f = 10000;

p = 1 * pi;

t = linspace(0, 10*pi, 100000);

s = [sin(2*pi*f*t); sin(2*pi*f*t + p)]'; %Setup signal

Ts = mean(diff(t));

Fs = 1/Ts;

Fn = Fs/2;

% fts = fft(s)/length(t); % Do transform

indxs = find(f>=9500 & f<=10500);

fts = goertzel(s, indxs);

Fv = linspace(0, 1, fix(length(t)/2)+1)*Fn;

Iv = 1:length(Fv);

[pk,ix] = max(abs(fts(Iv,:)));

phsmin = [min(angle(fts(ix(1),:))); max(angle(fts(ix(2),:)))];

phsdif = diff(phsmin); % Calculate Phase Difference

PhaseDiff = sprintf('Phase difference = %.3f rad', phsdif)