function PRC=prc(F,phases,x0,A) 
% plot phase-resetting curve (PRC) of system x'=F(t,x) + A delta(t) 
% at points given by the vector 'phases'. 
% 'x0' is a point on the limit cycle with zero phase 
% A is the strength of stimulation (row-vector) 
% use peaks of spikes to find the phase differences 
% Original code E.Izhikevich "Dynamical System in Neuroscience" chapter 10
% Modified by F. Piekniewski 2009 to work with Octave

% Example use:
% phases=0:0.1:2*pi;
% PRC=prc('F',phases,[-1 0],[0.3 0]);
%
% Where F:
% function dx = F(t,x); 
% z=x(1)+1i*x(2); 
% dz=(1+1i)*z-z*z*conj(z); 
% dx=[real(dz); imag(dz)];

T=phases(end); % is the period of the cycle 
tau = T/6000;	% time step of integration 
k=3;	% the number of cycles needed to determine the new phase 
PRC=[]; 
[tc,lc] = ode23s(F,0:tau:k*T,x0); % find limit cycle 
peak=1+find(lc(2:end-1,1)>lc(1:end-2,1)&lc(2:end-1,1)>=lc(3:end,1)); 
peak0 = tc(peak(end)); % the last peak is used for reference 
for i=1:length(phases)
   [m,j]=min(abs(phases(i)-tc)); 
   [t,x] = ode23s(F,phases(i):tau:k*T,lc(j,:)+A); % stimulate 
   peaks=1+find(x(2:end-1,1)>x(1:end-2,1)&x(2:end-1,1)>=x(3:end,1)); 
   PRC=[PRC, mod(T/2+peak0-t(peaks(end)),T)-T/2]; 
   subplot(2,1,2);drawnow; 
   plot(phases(1:length(PRC)),PRC); 
   xlabel('phase of stimulation');ylabel('phase difference'); 
   subplot(2,1,1); 
   plot(tc,lc(:,1),'r',t,x(:,1),t(peaks(end)),x(peaks(end),1),'ro'); 
   xlabel('time');ylabel('membrane potential');
end;

% This function substitutes ode23s which is not available in
% GNU Octave.
% coment for matlab (matlabs ode23s works MUCH faster)
%function [tc, lc]=ode23s(F,range,x0) % forward integration  
%tau=min(diff(range));
%lc(1,:)=x0;
%tc=range;
%range=range(2:end);
%for t=range 
%    lc(end+1,:)=lc(end,:)+tau*feval('F',t,lc(end,:)')';
%end;  