clear all
global G1 G2 V1 V2 r1 r2 i1 i2 E_Na E_K R L C Iin eps

%
% Parameter values for the Spike Burster IV-characteristics
%  1/r1<-G2
%
G1=1; G2=-0.5; V1=.5; V2=3;

r1=3; r2=-5; i1=.5; i2=.75;

E_Na = 0.5; E_K = 0.75;  


%
% Diagnostic plot
%
figure(1)
hold off
V=-1:.05:3;
i=-.5:.05:1.5;
plot(V,K(V+E_K),'r-',Na(i)+E_Na,i,'b-',[-1 3], [0 0], [0 0], [-1 1.5])


%
% Circuit simulation
%

R=0.01; L=3; C=.01; Iin=-6.5; eps=.001; % L=.5, 1.5, 2, 3 => 1, 2, 3, 4 spikes
% R=0.1; L=.5; C=.01; Iin=-7; eps=.001; % L=.3, .5 => 2, 3 spikes
%%%% R=0.1; L=.1; C=.01; Iin=-7; eps=.001; % close to Shilnikov orbit
% R--resistance, L--inductance, C--capacitance, 
% E--battery(Nernst potential)
% Increasing L to increase the number of spikes
% Increase Iin to suppress the spike bursts. 

%E_K=1;r2=-1;L=1;C=.001;eps=.0001;


T=1; % transcient time
% Numerical solutions
[t,x]=ode15s('SB',[0,8],[-0.2971;    6.3115;    0.1872]);
figure(2)
hold off
N=find(t>T); % Skip the first T interval solution 
plot(t(N:end),x(N:end,1),'r-',t(N:end),x(N:end,2),'b-',t(N:end),x(N:end,3),'g-')
hold on
legend('V_C','I_L','I_g')

figure(4)
hold off
subplot(3,1,1)
plot(t(N:end),x(N:end,1),'r-')
axis([1 8 -1 2.5])
ylabel('V_C')
subplot(3,1,2)
plot(t(N:end),x(N:end,3),'g-')
ylabel('I_g')
subplot(3,1,3)
plot(t(N:end),x(N:end,2),'b-')
ylabel('I_L')
xlabel('t')


figure(3)
hold off
plot3(x(N:end,1),x(N:end,2),x(N:end,3),'r-')
hold on


figure(1)
hold on
plot(x(N:end,1),x(N:end,3),'g-')

