clear all
global R L C E G1 G2 V1 V2  A T phi I_0 ex

G1=1; G2=-2; V1=.25; V2=.5;
% Nonlinear IV-characteristics: I = G1*V + (V>V1).*(V<=V_2).*(G2*(V-V1))+(V>V2).*(G2*(V2-V1));

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


A=0.5;T=10;phi=0;I_0=-1; 
% 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, T--period, ex--percentage of T that is the excitation interval

% Plot the nullclines 
%clear I V
%V=(-.5):.05:(2);
%figure(1)
%hold off
%plot(V,F(V),'r-')
%plot(V,-F(V+E)-A,'r-',V,-F(V+E)+A,'g-',V,V./R,'b-')
%axis([-.5 1 1.6 2.1])

% Numerical solutions
[t,x]=ode15s('FHN',[0,3*T],[-1;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-')

