GPregression_large_gap ;

%% Find all the point at a given frequency:
% Create new COCO problem - amplitude continuation
data = struct('x', x, 'Kiy', Kiy, 'par', par);
prob = coco_prob();

% Some options
prob = coco_set(prob, 'cont', 'h_min', 1e-4); % min step size
prob = coco_set(prob, 'cont', 'h_max', 0.05); % max step size
prob = coco_set(prob, 'cont', 'ItMX', 500); % max steps

% Add the test function
prob = coco_add_func(prob, 'evalGP', @eval_interpolant, data, 'zero', 'u0',[0.067,23.7,0.0089]) ;

% Tell COCO what the parameters are
prob = coco_add_pars(prob, '', 1:3, {'response_amp', 'forcing_freq', 'forcing_amp'});

% Add a slot function to monitor the progress of the continuation
% prob = coco_add_slot(prob, 'bifn_output', @bddat_slot_func, struct('ax', ax), 'bddat');

% Run the continuation
bd_amp = coco(prob, 'amp_sweep1', [], {'forcing_amp', 'response_amp'}, [1e-3, 0.3]);
xx = coco_bd_col(bd_amp, 'forcing_freq');
yy = coco_bd_col(bd_amp, 'forcing_amp');
zz = coco_bd_col(bd_amp, 'response_amp');

%% Do all frequency sweep using COCO
p0 = [zz',xx',yy'] ;
for i = 1:length(xx),
    
    try
        % Create new COCO problem - frequency continuation
        prob = coco_prob();
        
        % Some options
        prob = coco_set(prob, 'cont', 'h_min', 1e-4); % min step size
        prob = coco_set(prob, 'cont', 'h_max', 0.01); % max step size
        prob = coco_set(prob, 'cont', 'ItMX', 5000); % max steps
        
        % Add the test function
        prob = coco_add_func(prob, 'evalGP', @eval_interpolant, data, 'zero', 'u0', p0(i,:));
        
        % Tell COCO what the parameters are
        prob = coco_add_pars(prob, '', 1:3, {'response_amp', 'forcing_freq', 'forcing_amp'});

        % Run the continuation
        bd_freq = coco(prob, 'freq_sweep1', [], {'forcing_freq', 'response_amp'}, [23.65, 25.5]);
        
        % Plot the results
        xx = coco_bd_col(bd_freq, 'forcing_freq');
        yy = coco_bd_col(bd_freq, 'forcing_amp');
        zz = coco_bd_col(bd_freq, 'response_amp');
        plot(ax, xx, zz, 'linewidth',2.0, 'color', [128 128 128]/255); hold on;
        drawnow
    catch
        disp('Point missed') ;
    end
end