COCO

%% ========================================================================



% dydt = ...;  % define your continuation equation (vectorized, so don’t need to tell COCO otherwise)

x0 = zeros(22,1);
p0 = [0.0001];

p_range = [p0 10];   % parameter range for continuation

prob = coco_prob();        % define a new problem structure

ode_fcns = {@Coupled_Model, @Coupled_Modeldx,@Coupled_Modeldp};   % can append jacobian to this cell if you have it explicitly (ie. {dydt, dfdt, dfdp})

eq_path = {'path', 'to', 'isol2ep'};  % savdir/’run name’ (eg. /wd/data/path/to/isol2ep) by default sav dir is /wd/data/’run name’. input your path,

% helpful for saving different runs under different folder structures.


prob = ode_isol2ep(prob, '', ode_fcns{:}, x0, 'N(beta)', p0);  % initialize ‘initial sol to equilibrium point’ toolbox. x0, p0 from ode45/ode15s/etc…

% prob = coco_set(prob, 'cont', 'h_min', 0.1,'h0', 1,'h_max',10);    % stepsize preferences, can play around with these
%
% prob = coco_set(prob, 'cont', 'ItMX', [0 1e3], 'NPR', 1);      % max steps in run. print info every step.
%
% prob = coco_set(prob, 'corr', 'TOL', 1e-3,'ResTOL', 1e-3);        % correction tolerance, can play around with these



disp('  ')

disp('    initialising eq continuation.')

disp('-------------------------------------')

bd_eq = coco(prob, eq_path, [], 1, 'N(beta)', p_range);  % 1 is the solution dimension



%% ========================================================================

%///setup po continuation



labs = coco_bd_labs(bd_eq, 'HB');  % get the labels of the hopf bifurcations

bd_po = cell(1,length(labs));



if ~isempty(labs)

    disp('  ')

    disp('    initialising po continuation.')

    disp('-------------------------------------')

    disp('  ')

    for i = 1:length(labs) % you can parallelise continuation of multiple bifurcations to speed up the process if needed

        prob = coco_prob(); % new continuation problem structure

        prob = ode_HB2po(prob, '', eq_path, labs(i)); % HB2po toolbox, follows PO from HB

%        prob = coco_add_slot(prob, 'get_PO', @coco_savpo, [], 'bddat');  % needed to save complete PO data (but not min/max), can help you with @coco_savpo if needed

        prob = coco_add_event(prob, 'HB', 'BP', 'atlas.test.BP', 0); % stops continuation repeating itself if it rejoins the eq branch

        prob = coco_set(prob, 'corr', 'TOL', 1e-3,'ResTOL', 1e-3, 'ItMX', 15);

        prob = coco_set(prob, 'coll', 'NTST', 80);   % number of mesh points on which the collocation is performed

        prob = coco_set(prob, 'cont', 'NAdapt', 5, 'NPR', 1, 'PtMX', 1000);    % mesh adaption every 5 steps

        prob = coco_set(prob,...
            'cont', 'h_min', 0.1,'h0', 1,'h_max', 5,'ItMX', 1e3, 'al_max', 45);

        disp('done.')



        disp('  ')

        disp(['    starting po continuation at ', datestr(clock)])

        disp('--------------------------------------------------------')



        bd_po{i} = coco(prob,...
            {'path', 'to', strcat(sprintf('hb2po_%d', i))}, [], 1, {'N(beta)' 'po.period'}, p_range);  % output N(beta) and PO period during continuation

    end

else

    disp('no hopf points found for po continuation')

end



%% ========================================================================