% Plots the lower/upper values of limit cycle solution that can be superimposed on the
% equilibrium continuation curve

% We assume that limit cycle continuation has been called, e.g.,

%nmesh = 20;
%norder = 4;
%[x0,v0]=init_H_LC(@PlasmaMembraneOscill_ODEFunc, x1, p, [1], 1e-4, nmesh, norder);
%[xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);

% Then the following does the trick:

ODEDim = size(x,1)-1;
points = size(xlc, 2);
xx = xlc(1:end-2,:);
xx = reshape(xx, [ODEDim nmesh*norder+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]);