Кривая фазового отклика с MatCont (несколько кривых)

Question

Я не знаю, правильное ли здесь место, чтобы задать это, в любом случае, это вопрос программирования.

Я использую MatCont для расчета кривой фазового отклика модели Ходжкина-Хаксли, и я пытаюсь воспроизвести следующий рисунок.

из Кривые фазового отклика в нейробиологии, Натан У. Шультейс, с.261 .

В коде, во-первых, я нашел точку равновесия, по init_EP_EP точки бифуркации сохраняются в s. Я продолжаю кривую, где Хопф встречается в init_H_LC. Наконец, я продолжаю кривую с точки на предельном цикле (LC), используя функцию init_LC_LC. В каждом цикле я начинаю с одной и той же точки на LC и продолжаю больше точек (увеличивая MaxNumofPoints) и вычисляю кривую PRC.
Кажется, что-то не так, а также неэффективно.

global x v s h f opt

OPTIONS=[];
hls = HodgkinHuxley;
% find an equilibrium point
initial_condition = [-65.0 0.05293248525724958 0.3176769140606974 0.5961207535084603];
Istim=4.0;
[t,y] = ode45(hls{2},[0 1000],initial_condition,OPTIONS,Istim);
x0 = y(end,:)';

opt = contset;
opt = contset(opt, 'Singularities', 1);
opt = contset(opt, 'MaxStepsize', 1);
opt = contset(opt, 'MaxNumpoints', 1000);

ap = 1; %active parameter is Istim
Istim = 7.5;
[x1, v1] = init_EP_EP(@HodgkinHuxley, x0, Istim, ap);
[x,v,s,h,f] = cont(@equilibrium, x1,[], opt);

x1=x(1:4,s(3).index); % set initial state variables on H
Istim = x(end,s(3).index); % set the Inp where Hopf occures 
disp(Istim);

ap = 1;
ntst= 40; 
ncol= 4;
maxnumpoints=50; 

[x0,v0]=init_H_LC(@HodgkinHuxley,x1,Istim,ap,1e-6,ntst,ncol);
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-6);
opt = contset(opt,'MaxNumPoints',maxnumpoints);
[xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);

figure(2);
hold on;

for i1=1:4

    Istim = xlc(end,end);
    disp(Istim);

    [x0,v0]=init_LC_LC(@HodgkinHuxley,xlc,vlc,slc(end),Istim,ap,ntst,ncol);
    opt = contset(opt,'PRC',1);
    opt = contset(opt,'Adapt',1);
    opt = contset(opt,'MaxStepsize',5);
    opt = contset(opt,'FunTolerance',1e-6);
    opt = contset(opt,'VarTolerance',1e-6);
    opt = contset(opt,'Input',1);
    opt = contset(opt,'MaxNumPoints',40*i1);

    [xlc2,vlc2,slc2,hlc2,flc2]=cont(@limitcycle,x0,v0,opt);

    fvector=flc2(:,end);
    ind1 = ntst+2; ind2 =ntst+2+ncol*ntst; 
    PRC10=fvector(ind1:ind2);
    dPRC10=fvector((ind2+1):end);
    phi = linspace(0,1,length(PRC10));
    plot(phi, PRC10);
    legendInfo{i1} = ['Inp= ' num2str(xlc2(end,end))];

end

legend(legendInfo);

Файл HodgkinHuxley.m был загружен здесь на всякий случай.

Результат: