%here we estimate the slope of I_ion at the rest state

function HH_speed
global gnabar

dv = 0.0001;
V = [-12:0.0001:115];
gnabar = 120;  %nominal value

Ion_current = HH_rhs(V);
%find 3 zeros

ndx = find(Ion_current(1:end-1).*Ion_current(2:end)<=0);
   
 % slope at the rest point
 RM_inv = (Ion_current(ndx(1)+1)-Ion_current(ndx(1)))/dv
 
 RM = 1/RM_inv
Veq = -65;
figure(1)
plot(V,Ion_current,V(ndx(1)),0,'*',V(ndx(2)),0,'*',...
    V(ndx(3)),0,'*','linewidth',2);
xlabel('V','fontsize',20)
ylabel('-I_{ion}(V)','fontsize',20)
axis([-0.01 0.01 -0.004 0.004])


function Ion_current = HH_rhs(V)
global gnabar
gkbar = 36;

h = 0.5961;
n = 0.3177;

AM=.1*(25.-V)./(exp(.1*(25.-V))-1.);
BM=4.*exp(-V/18);
m = AM./(AM+BM);

gna=gnabar*m.^3.*h;
 
gk=gkbar*n.^4;

Vna = 115 ;
Vk =  -12 ;
Vl =  10.5988;
Icl=.3*(V-Vl);

Ion_current = -gna.*(V-Vna)-gk.*(V-Vk)-Icl;
