function HHdot=HHfunc(t,HHstate,flag,current)    
%flag
GNA = 120.0;	% mS/cm^2
GK = 36.0;	% mS/cm^2
GLEAK = 0.3;	% mS/cm^2

ENA = 115.0;	% mV
EK = -12.0;	% mV
ELEAK = 10.613;	% mV

C_M = 1.0;	% uF/cm^2

V = HHstate(1);     
M = HHstate(2);	% get values from last time step
H = HHstate(3);	% (hopefully this substitution makes 
N = HHstate(4);	% the following code a bit easier to read)

gNa = GNA * M^3 * H;
gK  = GK  * N^4;

mdot = alpha_m(V)*(1-M) - beta_m(V)*M;
hdot = alpha_h(V)*(1-H) - beta_h(V)*H;
ndot = alpha_n(V)*(1-N) - beta_n(V)*N;
vdot = (current*(t>5 & t<=30)...
   - gNa*(V-ENA) - gK*(V-EK) - GLEAK*(V-ELEAK))/C_M;  
HHdot=[vdot;mdot;hdot;ndot];
return

function rate = alpha_h(v)
rate = zeros(size(v));
rate = 0.07*exp(-v/20.0);
return

function rate = beta_h(v)
rate = zeros(size(v));
rate =  1.0 ./ (exp((-v+30.0)/10.0) + 1.0);
return

function rate = alpha_m(v)
rate = zeros(size(v));		% DEFAULT RATE TO ZERO     
idx = find(v~=25.0);			% HANDLE NORMAL (NON-SINGULAR) CASE
if(~isempty(idx))
   rate(idx) = 0.1*(-v(idx)+25.0) ./ (exp((-v(idx)+25.0)/10.0)-1.0);
end
idx = find(v == 25.0);		% HANDLE SIGULARITY AT V=25.0 
if(~isempty(idx))
   rate(idx) = 1.0/exp((-v(idx)+25.0)/10.0);
end
return

function rate =  beta_m(v)
  rate = zeros(size(v));
  rate = 4.0*exp(-v/18.0);
  return
  
function rate = alpha_n(v)
rate = zeros(size(v));
idx = find(v ==10.0);		% HANDLE SIGULARITY AT V=10.0 
if(~isempty(idx))
   rate(idx) = 0.1/exp((-v(idx)+10.0)/10.0);
end     
idx = find(v~=10.0);
if(~isempty(idx))
   rate(idx) = 0.01*(-v(idx)+10.0) ./ (exp((-v(idx)+10.0)/10.0) - 1.0);
end
return

function rate = beta_n(v)
  rate = zeros(size(v));
  rate = 0.125*exp(-v/80.0);
return