x1=x(1:2,s(5).index);p=[x(end,s(5).index);6];
[x0,v0]=init_H_LC(@MyODEFunc,x1,p,[1],1e-6,40,4);

opt = contset(opt,'MaxNumPoints', 100);
opt = contset(opt,'Multipliers',0);
opt = contset(opt,'Adapt',1);
opt = contset(opt,'MaxStepsize',5);
opt = contset(opt,'FunTolerance',1e-6);
opt = contset(opt,'VarTolerance',1e-3);

[xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);

figure(2)
clf
plotcycle(xlc,vlc,slc,[size(xlc,1) 1 2]);


figure(1)
hold on;

ODEDim = size(x,1)-1;
points = size(xlc, 2);
xx = xlc(1:end-2,:);
xx = reshape(xx, [ODEDim 40*4+1 points]);
Param = xlc(end,:);
xmax=squeeze(max(xx(1,:,:), [], 2));
xmin=squeeze(min(xx(1,:,:), [], 2));
cpl([Param; xmax'], vlc, slc, [1 2]);
cpl([Param; xmin'], vlc, slc, [1 2]);