clear

OPTIONS = [];

CaIinit = 0;
rhsfunc = PlasmaMembraneOscill_ODEFunc;
[t,y] = ode45(rhsfunc{2},[0 10],[0 0], OPTIONS, CaIinit);
x0 = y(end,:)';

figure(1)
subplot(311)
plot(t,y(:,1));
xlabel('[Ca^{2+}] = 0', 'FontSize', 18)
ylabel('Voltage, V(t)', 'FontSize', 18)

CaIinit = 0.2;
rhsfunc = PlasmaMembraneOscill_ODEFunc;
[t,y] = ode45(rhsfunc{2},[0 10],[0 0], OPTIONS, CaIinit);
subplot(312)
plot(t,y(:,1));
xlabel('[Ca^{2+}] = 0.2', 'FontSize', 18)
ylabel('Voltage, V(t)', 'FontSize', 18)


CaIinit = 0.4;
rhsfunc = PlasmaMembraneOscill_ODEFunc;
[t,y] = ode45(rhsfunc{2},[0 10],[0 0], OPTIONS, CaIinit);
subplot(313)
plot(t,y(:,1));
xlabel('[Ca^{2+}] = 0.4', 'FontSize', 18)
ylabel('Voltage, V(t)', 'FontSize', 18)


CaIinit = 0.0;
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxStepsize',1);
opt=contset(opt,'MaxNumpoints',100);


%----------- Do equilibrium continuation ------------------------------

[x1,v1]=init_EP_EP(@PlasmaMembraneOscill_ODEFunc, x0, [CaIinit], [1]);
[x,v,s,h,f]=cont(@equilibrium, x1,[], opt);

figure(2)
clf
subplot(211)
cpl(x,v,s,[(length(x0)+1) 1]);

open PlasmaMembraneBifDiag_withPeriod.fig

%----------- Do limit cycle continuation ------------------------------


singpt = 2;
x1=x(1:end-1,s(singpt).index);p=[x(end,s(singpt).index)];
nmesh = 100;
norder = 3;
[x0,v0]=init_H_LC(@PlasmaMembraneOscill_ODEFunc, x1, p, [1], 1e-4, nmesh, norder);

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

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

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


% ----------- Superimpose on equilibrium continuation ----------------

figure(2)
hold on;

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]);

axis([0 0.8 -90 40])


subplot(212)

xx = xlc(end-1,:);
Param = xlc(end,:);
cpl([Param; xx], vlc, slc, [1 2]);
cpl([Param; xx], vlc, slc, [1 2]);

axis([0 0.8 0 1.6])