function example_van_der_pol

% Example of computation of the Van der Pol oscillator'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-8,'AbsTol',1e-9); % set accuracy
init_cond=[1.0 0.0]; % initial condition

[T1 Y1]=ode45(@van_der_pol,[0 350],init_cond,ode_options); % integration over long time interval in order to achieve point on limit cycle
init_cond=Y1(end,:); % initial condition on the limit cycle

[period err]=periodas(@van_der_pol,7.0,init_cond,ode_options); % compute precise period of the limit cycle
[T Y PRC]=prc_fm(@van_der_pol,@jacobian_van_der_pol,'advanced',[period 2^10],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 = van_der_pol(t,y)

dydt(1,1)=y(2);
dydt(2,1)=1.5*y(2)*(1-(y(1))^2)-y(1);

function J = jacobian_van_der_pol(y)

J=[0 1; -3*y(1)*y(2)-1 1.5*(1-(y(1))^2)];