MatLab subroutine for computation of the phase response curve (PRC) of the stable limit cycle

Here you can download MatLab subroutines periodas.m, prc_fm.m and example_van_der_pol, example_hindmarsh_rose files. Put them in to your working directory and run one of the examples. The subroutines prepared for MatLab R2010b, but should work fine with older MatLab versions.

periodas.m

The subroutine periodas.m used to compute precise period of the limit cycle. The syntax is similar to the standard ode45 solver:

[tau err] = periodas(@sys_dynamic,appox_period,init_cond,options);

Here tau is a calculated period of the limit cycle, err is an error of the calculated period (should be same order as solver's accuracy).
@sys_dynamic is a function handle that evaluates the right side of the differential equations.
appox_period is an approximate value of the seeking period.
init_cond is initial conditions (should be on limit cycle).
options is a standard ode solver options which can be set via odeset function.

Since err usually can be neglected you can use such code: 

tau = periodas(@sys_dynamic,appox_period,init_cond,options);

prc_fm.m

The subroutine prc_fm.m realize algorithm of the PRC computation explained in the paper Computation of phase response curves via a direct method adapted to infinitesimal perturbations. The syntax is:

[T Y PRC] = prc_fm(@sys_dynamic,@jakob,method,period_points,init_cond,options);

Here T is a column vector of the time moments, Y is a matrix of the values of the limit cycle (each column corresponds to dynamical variable, while row corresponds to the time moments), PRC is a matrix of the values of the phase response curve (each column corresponds to the PRC component).
@sys_dynamic is a function handle that evaluates the right side of the differential equations.
@jacob is a function handle that evaluates the Jacobian matrix at the particular point on the limit cycle.
method must be equal to 'simple' or 'advanced'. The methods implement different multiplication procedure of auxiliary matrices. Both methods give the same accuracy of the computed PRC. The only difference is that 'simple' works with arbitrary number of intervals defined by int_number and it works slowly, while 'advanced' works faster but the number of intervals must be equal to \(2^m\) with \(m\) - natural number. If int_number is not of the form \(2^m\), when for the case of method = 'advanced' it will automatically decreased to the nearest \(2^m\). The computation time difference of two methods are most pronounced when int_number=2^12 or higher.
period_points = [tau int_number] is a two component vector where first term is precise period of the limit cycle and second term is the number of intervals in which the period is divided (note that equidistantly distributed node points where PRC will be computed is equal to int_number+1).
options is a standard ode solver options which can be set via odeset function. Note that even if number of intervals is small the accuracy of the computed PRC at the node points is equal to solver accuracy.

Mathematica notebook for computation of analytical expressions of Hamiltonian in the presence of high frequency expansion

The script toda_script.nb realizes an algorithm to compute the high-frequency expansion terms of the Hamiltonian. The algorithm is explained in the paper Flow-equation approach to quantum systems driven by an amplitude-modulated time-periodic force.