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;
%eps = 3;
%
% R--resistance, E--battery 
% eps = C/L with L--inductance, C--capacitance, 
% 


% Plot IV-characteristics
clear I V
V=0:.05:1.3;
I=F(V);
figure(1)
hold off
plot(V,I,'g--')

% Plot the nullclines 
figure(2)
hold off
plot(V-E,-F(V)-I_in,'r-', R*[-F(V(end))-I_in -I_in],[-F(V(end))-I_in -I_in],'b-')
%legend('V_C-nullcline','I_L-nullcline')

% Numerical solutions
[t,x]=ode15s('FHN',[0,5],[0;0]);


figure(3)
hold off
T=2;
N=find(t>T); % Skip the first T time
plot(t(N:end),x(N:end,1),'r-',t(N:end),x(N:end,2),'b-')
%plot(t,x(:,1),'r-',t,x(:,2),'b-')
legend('\itV_C','\itI_L','fontsize', 18)
xlabel('Time \it{t}','fontsize', 18)


hold on

% Orbit plot on the phase plane
figure(2)
hold on
plot(x(N:end,1),x(N:end,2),'mo')
%plot(x(:,1),x(:,2),'m-')
xlabel('\it{V_C}','fontsize', 18)
ylabel('\it{I_L}','fontsize', 18)

