%%%%%%%%%%%%%%% functions of HH Model %%%%%%%%%%%%%%%%%%
function xdot = hh_syst(t,x)

I = 0;     % applied current

v=x(1);
n=x(2);
m=x(3);
h=x(4);

an = 0.01*(-v+10)/(exp(-0.1*v+1)-1);
bn = 0.125*exp(-v/80);

am = 0.1*(-v+25)/(exp(-0.1*v+2.5)-1);
bm = 4*exp(-v/18);

ah = 0.07*exp(-v/20);
bh = 1/(exp(-0.1*v+3)+1);

dv = I + 36*n^4*(-12-v) ...
   + 120*m^3*h*(120-v) ...
   + 0.3*(10.6-v);

dn = an*(1-n)-bn*n;
dm = am*(1-m)-bm*m;
dh = ah*(1-h)-bh*h;

xdot = [dv dn dm dh]'; % column vector: '