function firefly
global Omega omega A
Omega=2; % stimulus frequency
omega=1; % firefly frequency
A=0.992;

y0(1)=0;
y0(2)=5*pi/4;

t=0;
h=0.06;

tint=100;
N=tint/h;

t_prev=0;
phase_diff_prev=0;

figure(1)
set(gcf,'DoubleBuffer','on');

clf;

subplot(2,1,1);
set(gca,'Xlim',[0.0,100],'Ylim',[-10,10]);
title('Phase difference (stimulus-firefly)');
subplot(2,1,2);
axis equal
set(gca,'Xlim',[-1.2,1.2],'Ylim',[-1.2,1.2]);
title('Red=Firefly phase; Blue=Stimulus phase')
prevPointStimulus = [];
prevPointFirefly = [];
hold on;
for i = 1:N
    tspan=[t t+h/2 t+h];
    [t1,y]=ode15s(@f1,tspan,y0);
    subplot(2,1,2)
    if ~isempty(prevPointFirefly);
      delete(prevPointFirefly);
    end
    
if ~isempty(prevPointStimulus);
      delete(prevPointStimulus);
    end
   
    
    prevPointStimulus = plot(cos(y(3,1)),sin(y(3,1)),'bo');
    prevPointFirefly = plot(cos(y(3,2)),sin(y(3,2)),'ro');
    drawnow;
subplot(2,1,1);

t=t1(3);
phase_diff=y0(1)-y0(2)

if (t_prev>0)
    line( [t_prev t],[phase_diff_prev,phase_diff]);
end


y0(1)=y(3,1);
y0(2)=y(3,2);

t_prev=t;
phase_diff_prev=phase_diff;

drawnow;
end
hold off


% -----------------------------------------------------------------
function dydt = f1(t,y)
global Omega omega A
dydt = [     Omega ;
             omega + A*sin(y(1)-y(2)) ]; 
% ----------------------------------------