clear all
global C E I_in r1 r2 i1 i2  eps A T phi ex

r1=1; r2=-2; i1=.5; i2=.75;
% Nonlinear 'S' IV-characteristics: v = r1*i +
% (i>i1).*(i<=i2).*(r2*(i-i1))+(i>i2).*(r2*(i2-i1)); 
% r1<-r2; i1<i2;

C=.1; E=.2;
% C--capacitance, E--battery


eps = .001; A=-1.25;T=1.5;phi=0;
% I = i1*A*sin(2*pi*t./T+phi), |i1*A|--amplitude, T--period, phi--phase

ex=50;
% For square wave forcing, T--period, ex--percentage of T that is the excitation interval

% Plot the nullclines 
clear i v 
i = (0):.05:(1.25);
v = S(i);
figure(1)
hold off
%plot(v,i,'r-')
plot(E+v,i,'r-',[-.5 1.25],[-i1*A, -i1*A],'b-')
%axis([-.5 1 1.6 2.1])

% Numerical solutions
[t,x]=ode15s('Scircuit',[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-')

