function example_hindmarsh_rose

% Example of computation of the Hindmarsh-rose neuron's phase response
% curve (PRC). The algorithm uses aditional subroutine files 'periodas.m'
% and 'prc_fm.m', which must be in your working directory.

ode_options = odeset('RelTol',1e-9,'AbsTol',1e-10); % set accuracy
init_cond=[-1.655538656276766 -12.712981825016282   1.346379929078901]; % initial condition on limit cycle

[period err]=periodas(@hindmarsh_rose,600,init_cond,ode_options); % compute precise period of the limit cycle
[T Y PRC]=prc_fm(@hindmarsh_rose,@jacobian_hindmarsh_rose,'advanced',[period 2^12],init_cond,ode_options); % compute the PRC

figure;
subplot(2,1,1);
plot(T(:),Y(:,1),'-b','LineWidth',1);
subplot(2,1,2);
plot(T(:),PRC(:,1),'-b','LineWidth',1);



function dydt = hindmarsh_rose(t,y)

dydt(1,1)=y(2)+3*y(1)^2-y(1)^3-y(3)+1.3;
dydt(2,1)=1-5*y(1)^2-y(2);
dydt(3,1)=0.001*(4*(y(1)+1.6)-y(3));

function J = jacobian_hindmarsh_rose(y)

J=[6*y(1)-3*y(1)^2 1 -1; -10*y(1) -1 0; 0.001*4 0 -0.001];