coco_use_recipes_toolbox bvp_v1 coll_v1 % Add bvp_v1 and coll_v1 toolboxes to search path

% The continuation problem structure encoded below corresponds to 151 zero
% functions (120 collocation conditions, 27 continuity conditions, and four
% boundary conditions) in terms of 154 continuation variables (150
% basepoint values, one interval length, and three problem parameters), a
% family of three monitor functions that evaluate to the problem
% parameters, and the corresponding inactive continuation parameters 'om',
% 'ro', and 'eps'. Its dimensional deficit equals 0. For 'eps'=0, there
% exists a periodic orbit for the reduced system, for which the rotation
% number is approximately three. A two-dimensional family of such 1:3
% phase-locked periodic orbits is obtained by releasing 'ro' and 'eps' and
% allowing these to vary during continuation.

p0 = [3.5; 0.35; 0];
[t x0] = ode45(@(t,x) lang(x, p0), [0 5.3], [0.3; 0; 0.4]); % Approximate periodic orbit

prob = coco_prob();
prob = bvp_isol2seg(prob, '', @lang, @lang_DFDX, @lang_DFDP,  ...
  t, x0, {'om' 'ro' 'eps'}, p0, @po_bc, @po_bc_DFDX);
prob = coco_set(prob, 'cont', 'h', 0.25, 'h_max', .5, 'MaxRes', 10);
prob = coco_set(prob, 'cont', 'almax', 35);
prob = coco_set(prob, 'cont', 'PtMX', 500, 'NPR', 50);
coco(prob, 'run1', [], 2, {'ro' 'eps'}, {[-2 2] [-0.5 0.5]});


figure(1)
clf
hold on
atlas = coco_bd_read('run1', 'atlas');
[fdata uidx] = coco_get_func_data(prob, 'bvp.seg.coll', 'data', 'uidx');
x0idx = uidx(fdata.x0_idx);
p_idx = uidx(fdata.p_idx);
plot_atlas(atlas,p_idx(2),p_idx(3),x0idx(1));
axis equal
axis tight
view(60,30)
drawnow

coco_use_recipes_toolbox % Remove bvp_v1 and coll_v1 toolboxes from search path