clear;

omega0 = 70*2*pi;
Omega = 1.5*omega0;
D = 1/2/300;
G = 0.008;
nl = 0.002665;

p0 = [omega0; Omega; D; G; nl];

%% EP Toolbox
prob = coco_prob();
prob = coco_set(prob, 'ode');
funcs = {@param_autonomous};
coll_args = [funcs, {[0; 0; 0], {'omega0' 'Omega' 'D' 'G' 'nl'}, p0}];
prob = ode_isol2ep(prob, '', coll_args{:});
cont_args = {1, {'Omega'}, {[1 1000]}};
bd1  = coco(prob, 'run', [], cont_args{:});

function y = param_autonomous(x,p)
    omega0 = p(1,:);
    Omega = p(2,:);
    D = p(3,:);
    G = p(4,:);
    nl = p(5,:);

    y(1,:) = x(2,:);
    y(2,:) = -2.*D.*x(2,:) - x(1,:).*(1 + (G - nl.*x(1,:).^2).*cos(Omega./omega0.*x(3,:)));
    y(3,:) = 1;
end

function y = param_autonomous2(x,p)
    omega0 = p(1,:);
    Omega = p(2,:);
    D = p(3,:);
    G = p(4,:);
    nl = p(5,:);

    y(1,:) = x(2,:);
    y(2,:) = -2.*D.*x(2,:) - x(1,:).*(1 + (G - nl.*x(1,:).^2).*cos(x(3,:)));
    y(3,:) = Omega./omega0;
end