clear all
global a b c eps R E I_in

a=2; b=-2; c=1.1;

%
% Nonlinear IV-characteristics: d*V.*(a*V.^2+b*V+c);
%

R=0.5; E=.4; I_in=-.6; % I_in=0;
eps = .001;
%
% R--resistance, E--battery 
% eps = C/L with L--inductance, C--capacitance, 
% 


% Numerical solutions
xinit=[0.8243; 0.2325];
epsa=.001;
epsb=1;
K=40;
epsinc=(epsb-epsa)/K;
Vextrema=[];
%options=odeset('e','on','RelTol',1e-12,'AbsTol',[1e-12;1e-12],'BDF', 'on');
options=odeset('e','on','BDF', 'on');
for k=0:K
   eps = epsa+k*epsinc;
   [t,x]=ode15s('FHN2',[0,5],xinit,odeset('BDF', 'on'));
   xinit=x(end,:);
   [t,x,aa,bb,cc]=ode15s('FHN2',[0,10],xinit,options);
   xinit=x(end,:);
   Vextrema=[Vextrema;[min([bb(end,1),bb(end-1,1)]),max([bb(end,1),bb(end-1,1)])]];
end
epsbif=epsa:epsinc:epsb;
plot(epsbif',Vextrema(:,1),epsbif',Vextrema(:,2))
hold on; 
plot([1 1.5], Vextrema(end,:))




