clear all
global R L C E a b c d  A T phi I_0 ex

a=4;b=1;c=2;d=0;
% Nonlinear IV-characteristics: I = a*V.*(V-b).*(V-c)+d;
% with d=0, V=0,b,c are the V-intersepts. 

R=0.1; L=.5; C=.1; E=.5;
% R--resistance, L--inductance, C--capacitance, E--battery
% Decreasing L to increase the number of spikes


A=2;T=20;phi=0;I_0=0; 
% I = A*sin(2*pi*t./T+phi)+I_0, A--amplitude, T--period, phi--phase
% shift, I_0--translation

ex=50;
% For square wave forcing, A--amplitude, T--period, ex--excitation interval

% Plot the nullclines 
%clear I V
%V=(-.5-E):.1:(2.5-E);
%figure(1)
%hold off
%plot(V,-F(V+E)-A,'r-',V,-F(V+E)+A,'g-',V,V./R,'b-')
%legend('V_C-nullcline, I_{in}=max','V_C-nullcline, I_{in}=min','I_L-nullcline')

% Numerical solutions
[t,x]=ode15s('FHN',[0,3*T],[0;0]);
figure(2)
hold off
N=find(t>T); % Skip the first period of forcing
plot(t(N:end),x(N:end,1),'r-',t(N:end),x(N:end,2),'b-')
hold on
%plot(t(N:end),Iin(t(N:end)),'g-') % use it when using Iin input
%plot(t(N:end),sqw(t(N:end)),'g-') % use it when using sqw input


%figure(1)
%hold on
%plot(x(N:end,1),x(N:end,2),'r-')

