hh.m

function V = hh(t,I,b)
x = [0 0 1];
g = [120 36 0.3];
V = zeros(1,length(t)); 
V(1)=50;
dt = 1e-3;
for i=2:length(t)
    a = exp(-(V(i-1)+55)/10)-1;
    if a == 0
        dn = (1-x(1)) * 0.1 - x(1) * (0.125*exp(-(V(i-1)+65)/80));
    else
        dn = (1-x(1)) * (-0.01*(V(i-1)+55)/a) - x(1) * (0.125*exp(-(V(i-1)+65)/80));
    end

    a = exp(-(V(i-1)+40)/10)-1;
    if a == 0
        dm = (1-x(2)) - x(2) * (4*exp(-(V(i-1)+65)/18));
    else
        dm = (1-x(2)) * (-0.1*(V(i-1)+40)/a) - x(2) * (4*exp(-(V(i-1)+65)/18));
    end
    dh = (1-x(3)) * (0.07*exp(-(V(i-1)+65)/20)) - x(3) / (exp(-(V(i-1)+35)/10)+1);
    dV =  - g(1)*x(2).^3*x(3)*(V(i-1)-50) - g(2)*x(1).^4*(V(i-1)+77) - g(3)*(V(i-1)+54.387)+b+I(i);
    x = x + dt*[dn dm dh];
    V(i)=V(i-1)+dt*dV;
end
end